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FOREWORD 


The  computer  program  SHIPFIT  was  originally  published  in  an  internal  NSWC 
Technical  Note  dated  23  April  1976.  Since  that  time  SHIPFIT  has  been  used  to 
determine  the  added  masses  of  ships  and  submersibles  in  a  number  of  ship 
response  studies  conducted  by  the  Explosion  Damage  Branch  (R14)  of  this  Center. 
Recent  requirements  for  horizontal  added  masses  and  the  desire  to  reference  the 
computer  program  in  publications  distributed  externally  have  brought  about  this 
updating  of  the  program  and  method. 

This  revision  was  sponsored  by  NAVSEA  (63R)  through  the  Undersea  Weapons, 
Warheads,  and  Fuzes  Block  program. 


Approved  by: 


KURT  F.  MUELLER,  Acting  Head 
Energetic  Materials  Division 
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CHAPTER  l 
INTRODUCTION 

The  approach  taken  in  this  report  follows  the  conventional  strip-theory 
technique  for  estimating  the  distribution  of  added  mass  along  the  length  of  a 
ship  undergoing  transverse  vibrational  motions.  The  flow  is  considered  to  be 
inviscid,  incompressible,  and  irrotational ,  and,  within  cross-sectional  planes 
spaced  length-wise  along  the  hull,  to  be  two-dimensional.  The  portion  of  the 
ship  between  two  planes  is  referred  to  as  a  strip  or  segment,  and  the 
intersection  of  a  plane  with  the  hull  is  called  a  section.  In  strip  theory  ship 
sections  are  considered  to  be  unvarying  in  form  within  each  segment. 

The  added  mass  of  a  ship  segment  in  this  two-dimensional  flow  field  is 
usually  a  close  approximation  to  the  added  mass  under  more  realistic  conditions 
of  flow  and  can  be  calculated  using  complex  variables  theory  and  conformal 
mapping  techniques.  Once  obtained,  this  approximation  may  be  improved  by 
multiplying  it  by  a  reduction  factor  equal  to  the  ratio  of  the  approximate  two- 
dimensional  and  exact  three-dimensional  added  masses  of  a  vibrating  ellipsoid  of 
revolution  fitted  to  the  ship  for  which  an  exact  solution  is  known.  The  added 
masses  calculated  by  SHIPFIT  do  not  contain  the  reduction  factor. 

The  theory  employed  in  SHIPFIT  is  adapted  from  the  1957  paper  of  Landweber 
and  Macagno,  in  which  formulas  were  developed  for  the  added  masses  of 
floating  and  submerged  bodies,  of  arbitrary  cross  section,  oscillating  in  both 
horizontal  and  vertical  directions.  The  added  masses  are  expressed  in  terms  of 
the  coefficients  of  a  conformal  transformation  that  maps  the  ship  cross- 
sectional  form  onto  a  circle  of  fixed  radius.  In  SHIPFIT,  the  forms  of  these 
expressions  are  slightly  changed  to  facilitate  their  coding.  For  totally 
submerged  bodies,  the  theory  is  somewhat  modified  so  that  the  cross  sections  of 
submersibles  can  be  viewed  in  their  usual  upright  orientation  and  to  facilitate 
the  evaluation  of  the  mapping  function  coefficients. 
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The  method  by  which  the  coefficients  of  the  conformal  transformation  are 
obtained  differs  from  previously  published  techniques.  In  SHIPFIT,  a  cubic 
spline  function  is  initially  fitted  to  the  data  points  describing  the  section 
form.  The  mapping  function  is  then  fitted  to  the  spline  function  by  means  of  a 
least  squares  Newton-Raphson  technique.  Benefits  of  the  method  appear  to 
include  a  considerable  improvement  in  both  speed  and  accuracy  over  other 
published  methods. 

In  addition  to  calculating  the  added  section  masses  per  unit  length, 

SHIPFIT  includes,  as  an  option,  a  provision  for  determining  the  equilibrium 
orientation  of  the  ship  with  or  without  user-specified  conditions  of  compartment 
flooding.  This  is  done  by  the  subroutine  ORIENT  prior  to  the  added  mass 
calculations . 


♦This  paragraph  refers  to  the  state-of-the-art  in  1976  when  SHIPFIT  was  first 
developed  (see  the  Foreword). 
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CHAPTER  2 
THEORY  AND  METHOD 


The  theory  applied  to  the  case  of  a  ship  (or  submersible)  floating  on  the 
surface  is  identical  to  that  of  Landweber  and  Macagno^  with  the  exception  of 
some  changes  in  notation.  Therefore,  only  results  and  a  brief  description  for 
this  case  will  be  included.  The  theory  employed  in  SHIPFIT  for  a  submerged 
vessel,  however,  differs  somewhat  from  that  of  Reference  1  and,  consequently, 
will  be  developed  more  fully. 

The  use  of  conformal  mapping  to  obtain  the  added  masses  of  objects  in 
two-dimensional  potential  flows  is  a  common  topic  in  hydrodynamics  texts. 
Rasically  it  consists  of  finding  a  transformation  that  will  map  the  streamlines 
of  the  object  onto  those  of  a  circle  or  flat  plane  where  the  complex  potential 
and  kinetic  energy  integral  are  more  easily  derived.  The  added  mass  is  then 
readily  obtained. 

In  Reference  l,  floating  and  submerged  hulls  are  handled  in  very  similar 

ways.  The  ship  on  the  surface  is  actually  treated  as  if  it  were  a  completely 

submerged  "double  hull"  formed  by  reflecting  the  submerged  portion  of  the  hull 
in  the  surface.  The  two  cases  of  interest,  then,  are  illustrated  in  Figure  1. 
The  two  halves  of  the  "double  hull"  are  imagined  to  vibrate  in  such  a  way  that 
the  free  surface  boundary  condition  is  satisfied  along  the  reflection  line 
exterior  to  the  ship  form.  The  appropriate  boundary  condition  for  most  cases  of 
interest  is  that  the  velocity  potential  vanish  along  the  free  surface  (this  is 
actually  the  high  frequency  condition).  For  vertical  motion  this  boundary 
condition  is  satisfied  if  both  halves  of  the  double  hull  vibrate  together  as  a 
single  rigid  form.  For  horizontal  motion,  the  vibrations  of  the  two  halves  are 

equal  but  in  opposite  directions.  In  both  cases,  the  added  mass  of  the  double 

Sin 1 1  is  twice  that  of  the  desired  value. 
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(a)  SUBMERGED  VESSEL 

(UNILATERAL  SYMMETRY) 
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For  the  two  cases  illustrated  in  Figure  1,  Landweber  and  Macagno  map  the 
streamlines  of  the  (single  or  double)  ship  form  onto  those  exterior  to  a  circle 
of  radius  r^  by  using  the  transformation 


Z  = 


> 


(1) 


where  Z  is  a  point  in  the  complex  plane  containing  the  ship  section,  hereafter 
called  the  Z-plane,  and  c  =  re  is  a  point  in  the  5-plane  containing  the 
circle  of  radius  r^.  The  representation  of  the  complex  potential  <p  +  i\|>  in 
the  {-plane  is  assumed  to  be  of  the  form 


q>  + 


(2) 


where  q>  and  \|i  are,  of  course,  the  velocity  potential  and  stream  function. 

By  substitution  of  the  representations  of  <p  and  from  Equation  (2)  into  the 
kinetic  energy  integral  of  a  fluid  at  rest  at  infinity,  i.e., 


T 


4>d<i»  , 


(3) 


where  p  is  the  fluid  density  and  the  integration  is  over  ail  boundaries  of  the 
fluid,  Landweber  and  Macagno  show  that  the  kinetic  energy  and  added  mass  per 
unit  length  can  be  expressed  as  rather  simple  sums  involving  the  coefficients 

i  i 

b  ,  b^,  ....  Finally,  the  b'  coefficients  are  expressed  in  terms  of  the 
a'  coefficients,  which  map  the  ship  contour  onto  the  circle,  by  satisfying  the 
boundary  condition  that  requires  the  normal  velocity  of  the  fluid  to  equal  that 
of  the  ship  at  the  ship  contour. 
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The  transformation  function  given  as  Equation  (1)  is  quite  general  and  may 
be  used  to  represent  an  arbitrary  shape.  It  is  desirable,  however,  to  build  the 
symmetry  properties  exhibited  by  the  forms  of  Figure  1  into  the  mapping 
function.  It  is  easily  verified  that  the  following  relationships  must  be 
satisfied : 


UNILATERAL  SYMMETRY 
(about  vertical  axis) 


Z(-c)  =  -z(c), 


(4) 


BILATERAL  SYMMETRY  ZU)  *  Z(c)  .  (5) 

Here  the  overbar  denotes  the  complex  conjugate.  By  substitution  of  Equation  (1) 
into  (4),  one  finds  that  unilateral  symmetry  requires  that  the  coefficients  of 
odd  powers  of  c  be  purely  real  and  the  coefficients  of  even  powers  of  c  be 
purely  imaginary.  In  the  same  manner,  bilateral  symmetry  may  be  shown  to 
require  real  odd  coefficients  and  vanishing  even  coefficients. 


Added  Masses  of  a  Ship  Floating  on  the  Surface 

For  the  bilaterally  symmetric  case  it  is  convenient  to  express  the  mapping 
function  as 


7  -  V  3  — Zn 

Z  "  L  3n  C 

n=  1 


where,  as  stated  above,  t'  2  coefficients  are  now  real.  This  is  the  form 
actually  used  in  SHIPFIT.  Also  in  SHIPFIT,  the  ship  form  is  mapped  onto  the 
unit  circle  (r^l)  rather  than  onto  a  circle  of  arbitrary  radius.  The 
differences  between  Equation  (6)  and  Equation  (1)  (applied  to  the  case  of 
bilateral  symmetry)  are  purely  notational.  The  correspondence  between 
coefficients  is  easily  shown  to  be 
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al  = 


a  = 
n 


a2n-3  / 


2n-3 


>  1 


(7) 


For  the  bilaterally  symmetric  case,  the  expressions  for  added  mass  derived 
by  Landweber  and  Macagno  can  now  be  written  in  the  form  that  appears  in  SHIPF1T 
by  substituting  Equa t ions ( 7 ) .  The  added  mass  per  unit  length  for  vertical 
oscillations  becomes 


\  =  11(5  3l^ai  +  a2^  ”  pS  ’ 


(8) 


where 


S  =  f  Z!  <3~2n)  .  (9) 

n=  1 

S  is  shown  by  Landweber  and  Macagno  to  be  simply  the  submerged  area  of  the  ship 
cross  section  (half  the  area  of  the  "double  hull").  This  result  is  obtained  by 
integrating  1/2  j)  xdy  over  the  double  hull  contour  after  substituting  parametric 
equations  for  x  and  y  in  terms  of  9. 

In  the  horizontal  direction  the  result  given  by  Landweber  and  Macagno  may 
be  written  as 


AH 


IT  Z  Z  Yjk  aj  +  l  ak+l  ’ 

j  =  l  k=  1 


(10) 
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CHAPTER  3 

USE  OF  THE  PROGRAM 


Input  data  is  supplied  to  SHIPF1T  through  three  separate  files:  DATA1 , 
DATA2,  and  INPUT.  DATA1  and  DATA2  are  presently  prepared  using  an  electronic 
digitizer  from  drawings  of  the  ship  cross  sections  and  the  ship  weight 
histogram,  respectively.  The  weight  distribution  is  needed  only  for  surface 
ships  or  surfaced  subrr.ersibles  and  only  if  a  computed  equilibrium  position  is 
desired . 


DATA1  File 


File  DATA1  contains  the  x  and  y  coordinates  relative  to  the  data  origin  of 
points  lying  along  the  ship  forms  and  the  locations  of  the  calculation  origins. 
Ship  forms  are  assumed  to  be  symmetric  and  either  the  left  or  right  half  of  a 
particular  cross  section  may  be  used.  The  order  of  the  points  on  DATA1  is 
important ;  an  improper  ordering  will  stop  the  program  execution  and  cause  the 

message  "IMPROPER  DATA  FORMAT  OR  CALCULATION  ORIGIN  FOR  SECTION  ...  STOP"  to 

be  printed.  The  first  point  of  each  section  is  the  calculation  origin  (see 
definition,  page  14).  The  second  and  successive  points  up  to  a  maximum  of  64 
must  begin  in  the  vicinity  of  the  lowest  point  of  intersection  of  the  ship  hull 
with  the  vertical  symmetry  axis  and  continue  in  steps  along  the  hull  in  either  a 
positive  or  negative  sense  of  rotation  about  the  calculation  origin.  For 
surface  ships,  the  last  point  should  define  the  freeboard  of  the  section.  For 
submersibles ,  the  last  point  should  be  in  the  vicinity  of  the  uppermost  point  of 
intersection  of  the  hull  with  the  vertical  axis.  The  data  are  read  in  the  main 
subprogram  of  SHIPFIT  by  the  list  directed  read  statement 


READ  (1,*)  IDIG ,  XDIG ,  YDIG 
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Determination  of  the  Ship  Equilibrium  Position 

As  an  option,  SHIPFIT  determines  the  orientation  of  a  ship  floating  on  the 
surface.  This  is  done  by  satisfying  the  two  equations 


X>  i  *  X>j 

1  j 

X  ®i*Bi  =  X!  WjZWj  ’ 


(44) 


(45) 


where  is  the  lumped  buoyancy  force  of  the  ith  ship  segment  centered  at 
7.^^  and  Wj  is  the  lumped  weight  from  the  jth  column  of  the  ship  weight 
histogram  centered  at  zWj •  Both  distances  are  measured  from  the  forward 
perpendicular  (forward-most  end  of  the  designer's  waterline).  A  provision  for 
additional  weights,  due,  for  example,  to  flooding  is  included.  Ship  segments 
pertaining  to  missing  cross  section  data  are  considered  neutral.  Equations  (44) 
and  (45)  are  satisfied  iteratively  in  subroutine  ORIENT  by  adjusting  the 
equation  of  the  ship  waterline. 
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For  sections  with  bilateral  symmetry,  all  but  the  first  two  coefficients 
are  assigned  inital  values  of  zero.  Initial  values  for  a^  and  a^  are  then 
obtained  from  the  half  beam  h  at  the  waterline  and  draft  d  of  the  section. 
Using  Equations  (33),  (34),  (37),  and  (38)  and  setting  values  of  6  to  0  and 
-w/2,  we  obtain 


ax  =  (d  +  h ) / 2 


a2  =  (d  -  h)/2 


(42) 


The  initial  coefficients  for  sections  with  unilateral  symmetry  are  obtained 
from  the  form  limits  along  the  y  axis,  say  yg  and  yg,  and  the  enclosed  area  A. 

A  is  calculated  from  the  spline  function  fit  to  the  section.  All  other 
coefficients  are  zero.  The  values  used  are 


(y. 


31  = 


_y_i!  + _ h _ _ 

*(y!  -  y 


a 


2 


(43) 


a 


3 


ys)/2 


This  set  of  initial  values  actually  amounts  to  a  four  parameter  fit  to  the  form 
as  a/t  is  then  identically  zero.  The  coefficients  specify  an  ellipse  enclosing 
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are  sought.  R($^)  is  a  functional  representation  of  the  ship  hull  relative 
to  the  calculation  origin.  SHIPFIT  uses  a  cubic  spline  function  for  this 
purpose  because  of  its  desirable  smoothing  properties,  speed  of  calculation,  and 
the  continuity  of  its  first  and  second  derivatives. 

The  Newton-Raphson  technique  is  an  iteration  scheme  in  which  the 
dE/Sa^  are  expanded  up  to  first  order  in  Taylor  series  about  the  initial 
set  of  coefficients.  We  obtain 


3E 

3a 


(a) 


n 


3  E  /  o . 

IT  (i  >  * 

n 


z 

k=l 


3?E 


3a.  3a 
k  n 


(a°)  [ a^  -  a£]  ,  n=l,2,...N, 


(40) 


where  the  underbar  denotes  the  full  vector  of  coefficients.  The  right  hand 
sides  of  Equation  (40)  are  equated  to  zero  and  the  set  of  N  linear  equations  are 
solved  for  a  new  set  of  coefficients  {a^},  which  is  used  to  replace  the 
old  in  the  next  iteration  step. 

SHIPFIT  employs  a  technique  to  accelerate  this  process  that  is  found  in  the 
NSWC/WOL  Library  subroutine  LSQSUB.  To  appreciate  this  feature  let  us  first 
expand  the  second  derivative  appearing  in  Equation  (40). 


3a,  3a 
k  n 


2 


z 

1 


2  ~\ 

34  34  3  4 

_ m  _ m  +  _ m 

3a,  3a  m  3  a,  3a 
k  n  k  n 


(41) 


where  4  =  r  -  R($  ).  The  technique  consists  of  deleting  the  second 

m  cm  cm  ^  ° 

derivative  term  on  the  right  hand  side  and  approximating  3  E/ 3a^3an  by  the 
remaining  terms  (twice  the  sums  of  the  products  of  the  first  derivatives).  This 
reduces  the  computation  time  significantly  as  second  derivatives  of  4^  need  not 
be  calculated.  If  the  process  fails  to  produce  convergence  in  a  given  number  of 
iterations,  SHIPFIT  repeats  the  process  with  the  exact  expression  for  the  second 
derivative. 
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For  an  approximately  correct  initial  set  of  coefficients  {a°}  and 
unilateral  symmetry,  a  half  plane  in  the  t-plane  will  map  onto  the  same  half 
plane  in  the  Z-plane.  Similarly,  for  the  bilaterally  symmetric  transformation, 
quadrants  will  map  onto  quadrants.  Such  mappings  are  suggested  in  Figure  1  by 
the  primed  and  unprimed  sets  of  points  and  image  points  A,B,C,D  and  A'.B'jC'.D'. 

i  6m 

SHIPF1T  creates  an  array  of  points,  e  ,m  =  1,2,  ...  M,  evenly  spaced  on  the 
unit  circle  and  spanning  the  quadrant  or  the  half  plane  of 

interest.  Images  (x  ,  y  }  of  these  points,  then,  lie  more  or  less  near  the 

sm  sm  r 

ship  form  in  the  Z-plane  depending  upon  the  correctness  of  the  initial  choices 

of  the  coefficients.  The  problem  consists  of  finding  a  way  of  improving  the 

initial  guesses  for  the  coefficients  so  that  the  points  {x  ,  y  }  lie  on 

sm  sm 

or  very  close  to  the  ship  form.  These  points  will  also  span  the  quadrant  or 
half  plane  of  interest  in  the  Z-plane,  but  the  spacing  between  the  points  in  the 
Z-plane  will  in  general  not  be  uniform. 

It  should  be  noLed  that  by  fixing  the  points  in  the  C-plane,  Equations 

(33)  and  (34)  become  very  simple.  As  noted  above,  the  functions  c^fm.n)  and 

Cy(m,n)  become  fixed  constants,  and  for  a  reasonable  range  of  n  values,  they 

may  be  tabulated  prior  to  the  task  of  improving  the  coefficient  values.  In  the 

2 

paper  by  von  Kercaek  and  Tuck  ,  points  in  the  Z-plane  were  held  fixed  and 
their  corresponding  image  points  in  the  c-plane  were  left  free.  The  location 
of  the  c-plane  points,  then,  required  the  solution  by  iteration  of  M 
transcendental  equations  involving  the  trigonometric  functions.  Computation 
times  using  their  method  are  thus  considerably  greater  than  those  for  the 
present  method. 

Adjustment  of  the  coefficients  to  increase  the  closeness  of  the  fit  is 

accomplished  in  SHIPFIT  by  means  of  a  least  squares  technique.  The  points  {x  , 

sin 

y  }  are  first  transformed  to  a  polar  representation,  (r  .  <t>  ),  about  the 

sm  ’  cm  cm 

calculation  origin.  Then,  by  means  of  a  Newton-Raphson  technique,  the  set  of 
coefficients  satisfying  the  equations 


11 

9a 

n 


(r  -  R(*  ) ) 2 )  =  0 ,  n=l , 2 , 

cm  cm 


.N 


(39 
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x  =  \  c  (m.n)  a 

sm  /  x  n 


00 

=  \  c  (m.n)  a 

L.  y  «  • 


In  the  case  of  unilateral  symmetry,  the  quantities  c  and  c  are  identified 

x  y 

from  Equations  (24)  and  (25)  as 


c  (m,n)  =< 


cos  [ (2-n)6  ] ,  n  odd 
m 


-sin  [(2-n)0  ],  n  even 
I  m 


c  (m.n)  = 

y 


sin  [ (2-n) 6  ] ,  n  odd 
m 


cos  [ (2-n) 9  ] ,  n  even 
m 


In  the  bilaterally  symmetric  case,  the  c  and  c  quantities  are  more  simply 

x  y 

expressed  as 


c  (m,n)  =  cos  [(3-2n)0  ] 
x  m 


c  (m,n)  »  sin[(3-2n)0  ] 
y  m 


Hence,  for  a  fixed  set  of  points  on  the  unit  circle  in  the  c-plane 
corresponding  to  a  fixed  set  of  angles,  the  quantities  c^  and  c  in 
Equations  (33)  and  (34)  are  also  fixed. 
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FIGURE  2.  COORDINATE  SYSTEMS  USED  IN  SHIPFIT 
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Fitting  of  the  Mapping  Function  Coefficients 

An  understanding  of  the  program  and  the  method  used  for  obtaining  the 
coefficients  requires  the  definitions  of  the  three  coordinate  systems  employed. 
These  are  shown  in  Figures  2a  and  2b.  As  they  are  related  through  simple 
translations,  they  are  described  by  definitions  of  their  origins. 

The  Data  Origin;  This  is  the  reference  origin  for  the  (x,y) 
coordinates  of  input  points  lying  along  the  section  that  define  its  shape.  The 
data  origin  must  always  be  along  the  vertical  axis,  and  for  ships  floating  on 
the  surface,  it  must  be  the  same  for  all  sections.  A  convenient  location  is  the 
intersection  of  the  vertical  symmetry  axis  with  the  base  line. 

The  Ship  Origin:  This  is  the  origin  of  the  Z  plane.  For  surface 
ships,  it  is  taken  as  the  intersection  of  the  vertical  axis  with  the  water 
line.  For  submerged  vessels,  it  is  located  at  the  geometric  center  of  the 
form. 


The  Calculation  Origin:  This  is  a  somewhat  arbitrary  point  chosen  so 
that,  in  a  polar  representation  about  this  origin,  the  radius  of  the  ship  form 
can  be  described  by  a  single  valued  function  of  the  angle.  This  is  a  relatively 
minor  constraint  on  the  generality  and  utility  of  SHIPFIT  as  a  calculation 
origin  can  almost  always  be  found. 

In  the  text  that  follows,  the  subscripts  d,  c,  and  s  will  be  used  to  reference 
variables  to  a  particular  coordinate  system  (data,  £alculation,  or  ship). 

The  objective  of  the  fitting  process  is  to  find  the  set  of  coefficients 

^®m 

(a^l,  that  will  make  points  lying  on  the  unit  circle  in  the  c-plane,  say  (e  , 

m=l,2 .  M},  the  images  under  the  mapping  transformation  of  points  lying 

on  the  submerged  portion  of  the  ship  hull  in  the  Z-plane,  say  {Zsm,  m=l,2..., 

M}  (reference  is  to  the  ship  origin). 

For  either  a  unilaterally  or  bilaterally  symmetric  ship  form,  the  x  and  y 
coordinates  of  the  points  in  the  Z-plane  can  be  expressed  generally  as  linear 
combinations  of  the  unknown  parameters  (a^),  i*e.. 
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iUa.l0  ,  j>2  and  even 

J+2  -  (29) 

b .  = 

J 

Ua  .  . „  ,  i >3  and  odd  . 

j+2  '  J— 

These  produce  an  added  mass  per  unit  length  for  a  submerged  vessel  vibrating 
horizontally  of  the  form 


ah  =  2*P  a1^a1  "  a3^  “  Ps'  » 


(3o: 


where  S'  is  again  given  by  Equation  (28). 

Comparison  of  the  formulas  for  A^  and  A^  with  those  obtained  by 
Landweber  and  Macagno  shows  the  forms  of  the  expressions  to  be  identical,  even 
though  the  symmetry  of  the  cross  section  is  now  about  the  vertical  axis  rather 
than  the  horizontal  or  imaginary  axis.  It  should  be  remembered,  however,  that 
the  coefficients  of  the  mapping  function  differ  in  the  two  approaches. 

Added  mass  coefficients  for  the  case  of  a  submerged  vessel  are  defined  as 


C 


V 


(3i: 


and 


C 


H 


irp(H/2)2 


(32; 


where  h  is  the  half  beam  of  the  submerged  form  and  H  is  the  full  height  from  top 
to  bottom  of  the  section. 
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The  b  coefficients  for  vertical  oscillations  are  found  by  inserting 
Equations  (18)  and  (24),  with  r=l,  into  Equation  (22)  and  matching  the 
coefficients  of  the  two  power  series.  We  obtain 


b. 

J 


real,  b^  = 


V°j*2  ' 


-iV(a3  +  a^ )  , 


j>2  and  even 
,  j_>3  and  odd  . 


(26) 


Substituting  these  values  into  Equation  (20),  we  obtain  the  added  mass  per  unit 
length  of  a  submerged  vessel  vibrating  vertically. 


Ay  =  2irp  a^(a^  +  a3)  -  pS '  , 


(27' 


where 


S' 


"I 


( 2-n)  a2 
n 


(28 


is  the  cross-sectional  area  of  the  submersible.  The  area  result  follows  upon 
substitution  of  Equations  (24)  and  (25)  into  &  x(3y/39)d9. 


In  a  similar  manner,  b  coefficients  for  oscillations  in  the  horizontal 
direction  are  found  to  be 


b0  =  iUa2  »  bl  =  U(a3  “  al)  ’ 
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S,  -  Sf  Y  nib 
H  ,,2  \  n 


These  results  are  obtained  by  simply  equating  T  from  Equation  (19)  and 
1/2  AyV2  or  1/2  ^U2. 

The  b  coefficients  in  Equations  (20)  and  (21)  may  be  expressed  in  terms  of 
the  a  coefficients  of  Equation  (15)  by  satisfying  the  boundary  conditions  on  the 
ship  contour.  Landweber  and  Macagno  show  these  to  be 


4>  =  -Vx 


for  vertical  oscillations  (velocity  V)  and 


\|»  =  Uy 


for  horizontal  oscillations  (velocity  U).  Here,  x  and  y  are  the  real  and 
imaginary  parts  of  Z  on  the  ship  contour,  which  may  be  expressed  as 


_  1  Vo  ,  i(2-n)0  +  **i(2-n)< 

"  2  L  3n(e 


n  odd 


-k  1 


■*  I 


i(2-n)0  -i(2-n)9 


i(2-n)0  .  -i(2-n)0> 


+; i  Y,  an(ei(2-n)e  -  e‘i(2-")9)  . 

n  odd 
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i.9 

where  is  a  complex  constant.  Substituting  c=re  ,  we  may  express  the 
real  and  imaginary  parts  of  this  as 


1 

2 


00 

I 

n=0 


1  /,  -in8  .  r  in9s 
—  (be  +  b  e  ) 

n  n  n 

r 


(17) 


CO 


(18) 


It  is  easily  shown*  that  this  form  of  the  complex  potential  yields  the  same 
result  for  the  kinetic  energy  as  obtained  by  Landweber  and  Macagno,  i.e., 


T'f^In|bni2-  <19) 

n=  1 


This  is  obtained  by  inserting  Equations  (17)  and  (18)  into  Equation  (3). 

The  added  masses  per  unit  length  for  submerged  hulls  vibrating  vertically 
with  velocity  V  or  horizontally  with  velocity  U  are,  respectively, 


(20) 


Proof  is  as  follows.  Let  <p'  =  <p  +  and  =  «|»  +  i|Iq.  Then 
^><p'd<i'  =  ^ipd\| i  +  d  i|» 

But  jidty  =  0  since  the  integration  may  proceed  around  the  ship  form,  to 
infinity  along  a  streamline,  around  the  circle  at  infinity  and  back  to  the 
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S. 


Added  Masses  of  a  Submerged  Vessel 


I 

P. 

i*. 

r. 

»*. 

.\ 

i 


For  the  case  of  unilateral  symmetry  (i.e.,  a  submersible  operating  at  a  depth 
where  surface  effects  are  negligible),  it  is  convenient  to  use  a  mapping 
function  of  the  form 


z  ■  Z  *»  '2~n  *  1  z 

n= 1  n-2 

n  odd  n  even 

where  the  coefficients  are  all  real. 

Landweber  and  Macagno  in  two  respects.  First,  Equation  (1)  does  not  include  an 
additive  constant  while  Equation  (15)  does  (a^).  This  is  useful  in  the 
fitting  process.  Second,  for  real  coefficients,  Equation  (1)  has  unilateral 
symmetry  about  the  real  axis  rather  than  about  the  imaginary  axis.  Thus,  to 
apply  the  results  of  Landweber  and  Macagno  to  submersibles,  the  cross  sections 
must  be  rotated  90  degrees.  It  is  desirable  to  have  a  theory  where  such  a 
rotation  is  unnecessary. 


2-n 


(15) 


This  form  differs  from  that  used  by 


b' 


Because  velocities  are  invariant  under  changes  of  the  complex  potential  by 
an  additive  constant,  the  transformation  for  the  complex  potential,  Equation 
(2),  can  be  written  more  generally  as 


bl  b2 

♦  i*  «  b  +  —  +  +  . 

C  ? 


(lb) 


t/ ■S'w->Sv y>Sy-i-Z 


ft 

. 

. 
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where 


s  ,1, 


a .  3a.,  j  >2 
J  J 


and  the  coefficients 
expressions : 


'  jk* 


j,k  >  l  are  computed  by  the  following 


( 


r  k'j 

(_  (  2  j- 1 )  (  2k- 1 )  V-"  _ 1 _  . 

4(k-j)  T2(  j+m)-l]  [2(k-m)-ir  *  ■* 

m=  1 

Tj*  *  ■  < 


2j-l 


z 

m=l 


_ 1 _ 

[  2  (  j  — m  )  —  1  ]  ^ 


j=k  . 


Added  mass  coefficients  are  defined 
ship  form  divided  by  the  added  mass  of  a 
floating  on  the  surface,  these  are  given 


as  the  ratio  of  the  added  mass  of  Che 
similar,  circular  form.  For  a  ship 
as 


C 


V 


a 


irph"/2 


( 


and 


CH 


2  ~  » 
irpd  /  2 


( 


where  h  is  half  beam  of  the  cross  section  at  the  waterline  and  d  is  the  ship 
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where  IDIG  is  the  sequential  integer  count  of  points  pertaining  to  a  particular 
cross  section.  Each  time  data  for  a  new  cross  section  appears  in  the  file 
(i.e.,  a  calculation  origin)  the  value  of  IDIG  must  be  reset  to  one.  Data  for  a 
total  of  24  cross  sections  may  appear  on  DATA1.  Finally,  points  should  be  more 
closely  spaced  where  the  ship  form  makes  sharp  turns  or  has  a  small  radius  of 
curvature.  This  is  because  the  derivatives  of  the  spline  function  fitted  to  the 
form  are  continuous.  An  example  of  DATAl  input  for  a  surface  ship  appears  in 
Figure  3.  Note  that  the  data  origin  is  set  at  the  intersection  of  the 
centerline  and  baseline,  and  that  the  calculation  origin  is  arbitrarily  located, 
but  chosen  so  that  R($)  is  not  highly  varying.  Point  coordinates  for  each 
section  may  be  multiplied  by  an  arbitrary  scale  factor  (see  INPUT  file 
instructions) . 


DATA2  File 

File  DATA2  contains  the  ship  weight  distribution.  It  will  be  used  if  the 
ship  is  floating  on  the  surface  and  if  the  equilibrium  calculation  is  selected 
by  setting  EQCAL  to  TRUE  on  the  INPUT  file;  otherwise,  DATA2  is  not  required. 
Data  are  read  in  the  subprogram  ORIENT  by  the  list  directed  read  statement 


READ  (2,*)  IDIG,  ZDIG ,  WZDIG 


where  IDIG  is  a  point  counter  (for  convenience  only),  ZDIG  is  the  position 
relative  to  the  forward  perpendicular  of  the  center  of  a  weight  histogram 
column,  and  WZDIG  is  the  column  height,  giving  weight  per  unit  length.  The  ZDIG 
scale  is  defined  as  zero  at  the  forward  perpendic ular  (bow)  and  positive  in  the 
direction  of  the  after  perpendicular  (stern).  The  ordering  of  points  on  DATA2 
is  not  important.  Any  point  sequence  may  be  used,  but  all  weight  histogram 
columns  should  be  represented,  up  to  a  maximum  of  48  columns  (the  program  only 
uses  those  columns  that  overlap  ship  segments  included  on  the  DATAl  file). 
Arbitrary  scale  multipliers  may  be  read  in  from  the  INPUT  file.  An  example  of 
DATA2  input  appears  in  Figure  4.  (Note:  the  number  of  weight  histogram  columns 
does  not  need  to  equal  the  number  of  ship  segments.) 
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FIGURE  3.  EXAMPLE  OF  DATA1  FILE  LISTING  FOR  A  SURFACE  SHIP  CROSS  SECTION 


FIGURE  4.  EXAMPLE  OF  DATA2  FILE  FOR  A  PARTICULAR  WEIGHT  HISTOGRAM 
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INPUT  File 


Data  on  the  INPUT  file  control  the  operation  of  the  program  and  provide  the 
values  of  certain  necessary  quantities.  The  structure  of  the  data  on  INPUT  will 
vary  according  to  the  intentions  of  the  user.  A  number  of  INPUT  records  are 
conditional  upon  the  values  of  variables  read  previously,  as  indicated  in  the 
READ  statement  structure  appearing  in  Figure  5.  Listed  below  are  variable 
definitions  in  the  order  in  which  they  are  read.  Line  designations  match  those 
in  Figure  5.  Many  variables  are  provided  with  default  values  that  are  used  if 
the  corresponding  input  field  is  left  blank.  These  are  underlined  in  Figure  5. 

A  list  of  these  variables  along  with  their  default  values  appears  in  Figure  6. 
Although  physical  dimensions  are  stated,  any  standard  consistent  system  may  be 
used  (Caution:  see  SPVOL,  Line  b). 


INPUT  File  Variable  Definitions  (Refer  to  lines  in  Figure  5) 

Line  a  (Title)  : 

TITLE  is  an  arbitrary  character  array  that  is  printed  on  each  page  of 
output  and  each  plot. 


Line  b  (Main  Program  Control  Parameters): 


SUBMERG  is  TRUE  if  the  ship  is  submerged  and  FALSE  if  it  is  floating. 

NNDFALT  is  the  default  value  of  NN,  the  number  of  coefficients 
calculated.  For  the  Ith  cross  section,  if  [nNA(I)J  >  NNDFALT  (see 
NNA(I)  description,  Line  c)  the  first  NN  =  NNDFALT  coefficients  are 
found  in  a  single  iteration  process;  then  NN  is  increased  by  1  per 
iteration  process,  using  the  previous  set  as  initial  values,  until  NN 
=  | NNA ( I ) |  parameters  have  been  determined.  If  left  blank,  NNDFALT=8 
is  assumed. 

PLOT  =  TRUE  results  in  the  plotting  subroutine  SHIPLOT  being  called 
(see  Line  d).  Blank  is  interpreted  as  FALSE. 
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READ  STATEMENTS:* 

a.  READ  5,  (TITLE(J),  J  =  1,  8) 

b.  READ  6.  SUBMERG.  NNDFALT.  PLOT.  AIPRINT.  SPVOL  ** 

c.  READ  15,  <IMS(I).  FRAME(I).  SCALE(I).  NNA(I),  NDA(I).  I  =  1,  III) 

Note:  III  is  the  total  number  of  ship  cross  sections  appearing  on  DATA1  (lll<24) 

d.  If  PLOT  is  TRUE:  READ  20.  (IPLOT(I).  I  =  1. 10).  MULTI.  FMAG.  XEND.  YEND 

e.  If  the  ship  is  floating,  i.e.,  SUBMERG  is  FALSE: 

el.  READ  35,  LBP,  NSEGS,  EQCAL 

e2.  If  EQCAL  is  FALSE:  READ  40,  WATERLN 

e3.  If  EQCAL  is  TRUE:  READ  110,  COLWID,  WSCALE.  ZSCALE.  FLOODED 

e3i.  If  FLOODED  is  TRUE:  READ  135,  NXTRA,  (ZXTRA(J),  WXTRA(J),  J  *  1,  NXTRA) 

f.  READ  75,  NCASES-  (NSPICK(J).  J  =  1,  NCASES) 

g.  If  NNAOXO  for  the  Ith  section  being  processed:  READ  105,  IN,  (AI(J),  J  -  1,  IN) 

Note:  This  read  statement  is  executed  once  for  each  section  having  a  negative  NNA(  * )  value  read  in  through 
statement  c  above. 


FORMAT  STATEMENTS:* 


5  FORMAT  (8A10) 

6  FORMAT  (L10, 110,  2L10,  E10.0) 
15  FORMAT  (110,  2E10.0,  2110) 

20  FORMAT  (1011,  L10,  3E10.0) 

35  FORMAT  (E10.0, 110,  L10) 


40  FORMAT  (E10.0) 

75  FORMAT  (1615) 

105  FORMAT  (I10/(8E10.0) 
110  FORMAT  (3E10.0,  L10) 
135  FORMAT  (I10/(2E10.0)) 


•These  statements  indicate  the  logical  structure  of  INPUT  file  data  and  may  differ  from  statements  actually  appearing 
in  the  program. 

••Underlined  variables  have  default  values  that  are  used  if  the  data  field  is  left  blank  (see  definitions). 


FIGURE  5.  INPUT  FILE  READ  STATEMENT  STRUCTURE  AND  FORMATS 
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VARIABLE 

DEFAULT  VALUE 

CONSEQUENCE 

SUBMERG 

FALSE 

The  ship  is  assumed  to  be  on  the  surface  • 

NNDFALT 

8 

Eight  mapping  function  coefficients  are  calculated  - 

PLOT 

FALSE 

SHIPLOT  (plotting  subroutine)  is  not  called . 

AIPRINT 

FALSE 

Computed  coefficient  values  are  not  written  on  AIFILE  • 

SPVOL 

35  FT3/TON 

These  units  are  assumed.  Saltwater  is  assumed . 

SCALEO) 

1. 

DATA1  section  dimensions  are  not  rescaled. 

NNA(I) 

0 

NNDFALT  parameters  are  calculated. 

NDA(f) 

0 

Fast-fitting  scheme  for  coefficients  is  used. 

MULTI 

FALSE 

All  chosen  sections  appear  on  a  single  plot . 

FMAG 

1.0 

Y  axis  of  plot  is  six  inches  long. 

XEND 

30. 

X-axis  length  from  ship  origin  is  30  feet. 

YEND 

30. 

Y-axis  length  from  ship  orgin  is  30  feet. 

WSCALE 

1. 

Weight/length  data  on  DATA2  are  not  rescaled. 

Z  SCALE 

1. 

Distance  data  on  DATA2  are  not  rescaled. 

FLOODED 

FALSE 

Ship  is  watertight.  No  additional  weights  read. 

NCASES 

0 

Added  masses  are  calculated  for  all  sections. 

FIGURE  6.  DEFAULT  VALUES  FOR  SHIPFIT  INPUT  VARIABLES 
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AIPRINT  =  TRUE  results  in  output  printed  on  A1FILE  consisting  of  fitted 
parameter  values  (see  Line  f  below).  Blank  is  interpreted  as  FALSE. 

SPVOL  is  the  specific  volume  of  seawater,  usually  expressed  in 
i  3 

ft  /ton.  If  left  blank  (i.e.,  0),  a  default  value  of  35  ft  /ton, 

that  is  appropriate  for  standard  seawater,  is  used.  (Fresh  water  value 

is  36  f t^/ton. ) 


Line  c  (Control  Parameters  for  Specific  Sections): 


NS(I)  is  the  station  number  of  the  section.  Station  0  is  usually 
located  at  the  forward  perpendicular  (see  Line  f  below). 

FRAME(l)  is  the  frame  number  of  the  section.  Frame  numbers  are  not 
used  in  calculations,  but  appear  in  the  printout  and  plots.' 

SCALE ( I )  is  a  length  scale  factor  that  multiplies  both  DATA1  input 
coordinates.  Blank  implies  a  default  value  of  one. 

NNA(I)  is  the  number  of  mapping  function  coefficients  to  be  found  for 
the  section.  The  program  storage  limits  allow  for  a  maximum  of  24 
parameters.  If  NNA(I)  =  0,  NNDFALT  is  used.  If  NNA(I)  <  0,  initial 
values  for  the  coefficients  are  read  in  by  the  user  (see  Line  g  below). 

NDA(l)  =  0,  if  the  fast  approximation  of  the  Newton-Raphson  iteration 
process  described  above  is  to  be  used  initially.  If  convergence  does 
not  occur  in  50  iterations,  the  exact  method  follows  automatically. 

Any  other  value  will  result  in  the  exclusive  use  of  the  exact  method. 


Line  d  (Plotting  Control  Parameters): 

tPLOT(I)  is  assigned  an  integer  value  of  1  through  9  that  controls  the 
darkness  of  the  plotted  line  or  symbol  for  the  Ith  plot  category 
described  below.  A  zero  value  causes  the  Ith  plot  category  to  be 
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skipped.  Lower  values  produce  lighter  lines,  and  higher  values  produce 
darker  lines.  A  line  of  standard  thickness  is  produced  by  the  value  4. 
Available  options  are: 

IPLOT(l)  *  0  ...  Draw  fitted  spline  function  curve, 

IPL0T(2)  **  0  . . .  Plot  input  points  defining  section  (plotting  symbol 
is  a  circle) , 

IPL0T(3)  ^  0  ...  Draw  fitted  mapping  function  curve, 

IPL0T(4)  *  0  ...  Plot  Z-plane  images  of  points  on  unit  circle 
(plotting  symbol  is  a  triangle). 

(Variables  IPLOT(5)  through  IPLOT  (10)  are  reserved  for  future  use.) 

MULTI  *  TRUE  results  in  multiple  plots  —  one  for  each  section  picked 
in  Line  f.  FALSE  results  in  all  picked  sections  being  plotted  on  a 
single  plot  (see  example  in  Figure  13,  p.  39). 

FMAG  is  a  magnification  factor  that  can  be  used  to  shrink  or  increase 
the  size  of  the  plots,  e.g.  ,  FMAG  =  .5  produces  a  half  size  plot,  while 
the  plot  resulting  from  FMAG  “  1.5  would  be  50  percent  larger  than  the 
standard  size  which  has  a  y-axis  length  of  six  inches.  The  default 
value  is  1.0. 

XEND  is  the  distance  in  feet  from  the  ship  origin  to  the  end  of  the 
plotted  x  axis.  It  must  be  a  multiple  of  10.  The  default  value  is  30. 

YEND  is  the  distance  in  feet  from  the  skip  origin  to  the  end  of  the 
plotted  y  axis.  It  must  be  a  multiple  of  10.  The  default  value  is  30. 


Line  el  (Information  Required  If  the  Ship  Is  Floating  on  the  Surface): 

LBP  is  the  ship  length  between  perpendiculars  (length  at  the  designer's 
waterline)  in  feet.  LBP  is  REAL. 

NSEGS  is  the  total  number  of  ship  segments  (used  to  calculate  segment 
length,  LBP/NSEGS).  Default  value  is  20. 
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EQCAL  =  TRUE,  if  Che  equilibrium  position  is  Co  be  calculated.  Blank 
is  interpreted  as  FALSE. 


Line  e2  (Waterline  to  be  Used  for  Surfaced  Ship): 


WATERLN  is  the  waterline  position  relative  to  the  data  origin  in  feet. 


Line  e3  (Parameters  for  Determining  Weight  Histogram  from  DATA2  Input): 


C0LU1D  is  the  column  width  of  the  weight  histogram  in  feet. 

WSCALE  is  a  scale  factor  that  multiplies  the  weight  per  unit  length 
data  read  in  from  DATA2.  Blank  implies  a  default  value  of  one. 

Z SCALE  is  a  scale  factor  that  multiplies  the  weight  histogram  column 
position  data.  Blank  implies  a  default  value  of  one. 

FLOODED  is  TRUE  if  the  ship  has  been  flooded  and  extra  lumped  weights 
are  to  be  read  from  the  INPUT  file. 


Line  e3i  (Additional  Weight  Distribution): 

NXTRA  is  the  total  number  of  extra  lumped  weights  and  weight  positions 
to  be  read.  NXTRA  <  24. 

ZXTRA  is  an  array  of  extra  lumped  weight  positions  in  feet  relative  to 
the  forward  perpendicular  with  the  same  sign  convention  as  in  the 
weight  histogram  (zero  at  the  forward  perpendicular  and  positive  in  the 
direction  of  the  after  perpendicular). 

WXTRA  is  an  array  of  lumped  weights  in  (long)  tons. 
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Line  f  (Sections  Chosen  For  Added  Maas  Computations): 


NCASES  is  the  number  of  cross  sections  of  DATA1  to  be  processed.  If 
NCASES  =  0  (Blank),  all  forms  are  processed  in  the  order  of  their 
appearance  on  DATA1. 

NSPICK  is  an  array  of  section  station  numbers  in  the  order  in  which 
processing  is  to  occur. 


Line  g  (Optional  Input  of  Initial  Coefficient  Values): 

IN  is  the  number  of  initial  coefficient  values  to  be  read  in.  (If  the 
number  NN  of  coefficients  to  be  calculated  is  greater  than  IN,  the 
additional  NN-IN  coefficients  are  given  initial  values  of  zero.) 

AI  is  an  array  of  initial  coefficient  values.  Note  that  if  NNA(I) 

>_  0  for  the  Ith  section,  initial  coefficient  values  are  generated 
internally  by  the  program. 


INPUT  File  Example 


Figure  7  shows  the  INPUT  file  records  for  a  sample  problem  concerning 

a  class  of  destroyers.  (The  examples  shown  in  Figures  3  and  4  also 

pertain  to  this  ship  design.)  In  Figure  7,  considerable  use  is  made  of  the 

default  features  built  into  SHIPFIT.  The  blank  fields  on  Card  2  indicate  that 

the  ship  is  floating  (SUBMERG  is  FALSE),  the  default  value  for  the  number  of 

mapping  function  coefficients  is  8  (NNDFALT  =  8),  no  printing  of  final 

coefficient  values  on  AIFILE  is  desired  (AIPRINT  if  FALSE),  and  seawater  is 

3 

assumed  (SPVOL  =  35  ft  /ton).  Cards  3  through  17  give  station  and  frame 
numbers.  Here  SCALE(I),  NNA(I),  and  NDA(I)  take  their  default  values  indicated 
in  Figure  6.  Card  18  is  needed  because  PLOT  is  TRUE  in  Card  2.  Card  18 
indicates  that  spline  functions  will  be  drawn  with  a  standard  line  thickness 
(IPLOT(l)  =  4  in  column  1)  and  input  points  will  be  plotted  with  somewhat 
heavier  symbols  (IPL0T(2)  *  6  in  column  2).  Individual  plots  will  be  prepared 
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for  each  section  selected  (MULTI  is  TRUE),  and  a  75  percent  reduced  plot  will  be 
produced  with  x  and  y  axis  lengths  of  30  and  20  feet,  respectively.  Card  19 
shows  the  length  between  perpendiculars  LBP,  with  20  segments  assumed  by 
default.  Because  EQCAL  is  blank  (FALSE)  on  Card  19,  Card  20  must  supply  the 
position  of  the  waterline  relative  to  the  data  origin  used  for  the  DATA1  file. 
Finally,  Card  21  says  one  section  is  to  be  analyzed  (NCASES=1)  and  the  selected 
section  is  at  station  number  3  (NSPICK(l)  =  3). 

The  output  produced  by  this  sample  problem  is  shown  in  Figures  8,  9,  and  11 
and  is  discussed  below. 


Output 

SHIPFIT  produces  three  different  kinds  of  output:  printout,  plots 
(optional),  and  a  file  of  calculated  mapping  function  coefficient  values  called 
A1FILE  (optional).  The  principal  form  of  output  is  the  printout.  A  page  of 
information  detailing  the  various  aspects  of  the  added  mass  calculation  is 
printed  out  for  each  selected  section.  An  example  produced  by  the  input  of 
Figure  3,  4,  and  7  is  shown  in  Figure  8.  Results  included  are  the  residuals  of 
the  fit  of  the  spline  function  to  the  input  points,  the  RMS  error  of  the  spline 
fit,  initial  and  final  values  of  the  mapping  function  coeffients  showing 
individual  and  overall  quadratic  mean  errors  for  the  fitted  coefficients,  the 
calculated  ship  form  as  represented  by  the  Z-plane  images  of  points  equally 
spaced  about  the  unit  circle,  differences  of  colinear  radii  drawn  from  the 
calculation  origin  to  the  image  points  and  to  the  spline  function  curve,  the 
submerged  beam  and  area,  the  buoyancy  coefficient  (this  is  the  ratio  of  the 
submerged  area  to  the  area  within  a  semicircle  with  a  diameter  equal  to  the 
submerged  beam),  the  draft,  and  horizontal  and  vertical  added  masses  per  unit 
length  and  added  mass  coefficients.  The  distinction  between  the  "best"  and 
"check"  values  of  the  vertical  added  mass  is  explained  ’n  the  following  section 
on  program  operational  guidelines. 

If  the  ship  is  floating  on  the  surface,  the  pages  describing  the  added  mass 
calculations  are  preceded  by  a  page  of  ship  orientation  information.  Formats  of 
this  page  differ,  depending  upon  whether  the  waterline  is  read  in  from  the  INPUT 
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file  or  calculated  by  balancing  the  forces  and  torques  that  bear  on  the  ship. 
Figures  9  and  10  show  these  two  forms  of  the  output.  Figure  9  is  the  result  of 
the  input  given  in  Figure  7.  An  additional  column  of  extra  weights  is  included 
in  Figure  10  if  these  are  present  in  the  INPUT  file. 

Forms  of  the  plotted  output  are  shown  in  Figures  11,  12,  and  13.  Figures 
11  and  12  show  the  individual  plots  that  are  produced  when  the  variable  MULTI  is 
TRUE,  while  Figure  13  is  the  type  of  result  obtained  when  MULTI  is  FALSE. 

Figure  11  was  produced  by  the  INPUT  file  shown  in  Figure  7.  It  is  a  plot  of  the 
input  points  and  fitted  spline  function  for  the  cross  section  located  at  station 
3.  Figure  12  shows,  in  addition  to  the  spline  function,  plots  of  the  image 
points  (triangles)  and  the  fitted  mapping  function. 


Program  Operation  Guidelines 

It  is  apparent  from  the  description  of  the  INPUT  file  that  SHIPFIT  allows 
some  flexibility  in  the  manner  in  which  the  parameter  fitting  is  carried  out. 
Although  the  program-assigned  set  of  initial  coefficient  values  given  in 
Equations  (42)  and  (43)  usually  are  sufficient  to  produce  convergence,  an 
occasional  section  form  may  be  encountered  where  these  conditions  fail  to 
produce  convergence.  For  such  forms  it  may  then  be  necessary  to  set  NNDFALT  to 
a  value  smaller  than  the  desired  number  of  coefficients  NNA  (II)  where 
convergence  is  guaranteed  (such  as  6  or  8)  and  obtain  the  desired  NNA  (II) 
coefficients  incrementally.  A  single  form  that  is  difficult  to  fit  can  be 
worked  on  independently  by  using  the  NSPICK  array.  The  incremental  fitting  of 
parameters  can  be  done  either  automatically  as  described  in  the  NNDFALT 
definition  (Line  b) ,  or  in  a  single-step-at-a-time  procedure  using  the  AI  array, 
and  where  NNDFALT  is  incremented  in  successive  runs. 

A  problem  in  fitting  may  occur  along  concave  portions  of  a  particular  ship 
contour  as  shown  in  Figure  12.  Points  that  are  equally  spaced  on  the  unit 
circle  tend  to  map  into  points  that  are  more  tightly  spaced  along  convex  lengths 
of  the  ship  contour  and  less  closely  spaced  along  concave  parts.  This  presents 
a  fitting  problem  that  could  be  remedied  by  locating  more  points  in  regions  of 
the  unit  circle  that  map  into  the  concave  regions  of  the  contour.  Presently, 
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1  *0*9X*ADDED  MASS  /L  (BEST)*7X,F 1 2. 2, 5X,F 1 2. 2) 

PRINT  294, AV2 

294  FORMAT (*0*9X* ADDED  MASS  /L  (CHECK)*6X,F12.2) 

PRINT  295,CV,CH 

295  FORMAT (*0*9X*ADDED  MASS  COEFF 1 C I ENT*6X,F 1 2. 5, 5X,F 1 2. 5) 
I F (PLOT)  CALL  SHIPL0T(2) 

300  CONTINUE 

400  IF (PLOT)  CALL  SHIPL0T(3) 

STOP 

END 
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264  FORMAT (10X,2F7.2,F8.4,3X,2F7.2,F8.4,3X,2F7.2,F8.4,3X,2F7.2,F8.4) 

C 

C**»  CALCULATE  CROSS  SECTIONAL  AREAS  FROM  MAPPING  FUNCTION  COEFFICIENTS 
SUM=0. 

270  DO  275  K=1,NN 

ABK=ALPHA-BETA*FLOAT ( K ) 

SUM=SUM+A8K*AA(K, 1 )**2 
275  CONTINUE 

AREAMF=PI *SUM 

IF ( .NOT.SUBMERG)  AREAMF =AREAMF/2 . 

C 

C**»  CALCULATE  HYDRODYNAMIC  COEFFICIENTS  FOR  SECTION 
BEAM=BD2*2. 

CBUOY=2.*AREA( I  I )/PI/BD2**2 
IF(SUBMERG)  CBU0Y=CBU0Y/2. 

PRINT  280, BEAM, AREA ( I  I ),CBU0Y 

280  FORMAT ( *0*9X*SUBMERGED  BEAM  =*F7.2//10X*SUBMERGED  AREA  =*F7.2// 

1  10X*8UOYANCY  COEFFICIENT  =*F10.5) 

I F ( SUBMERG)  GO  TO  290 

AVI  =DENSI  TY*(PI  *AA(  1  , 1  )  *  <  A  A  <  1 , 1  )  +AA  (2,1  ) )  -AREA  ( ID) 

AV2=DENS ITY*(PI *AA( 1 , 1 )*(AA(1,1 ) +AA{ 2, 1) ) -AREAMF ) 

AVCYL=DENS I TY*PI *BD2*#2/2. 

CV=AV1 /AVCYL 
SUM1 =SUM2=SUM3=0. 

N1 =NN-1 

DO  286  I V=2,N1 

SUM1=SUM1+AA( IV+1,1 )»CC( 1 , IV) 

SUM2=SUM2+AA( IV+1 )*#2*CC( I  V, I  V) 

I F { IV.EQ.2)  GO  TO  286 
I  V 1  =  I  V  —  1 

DO  284  I U=2, I VI 

SUM?=SUM3+AA( IU+1 ,1 )*AA( IV+1,1 )*CC( IU, IV) 

284  CONTINUE 
286  CONTINUE 

AH=8. *DENS I TY/PI *( ( AA (2, 1 )-AA (1,1) )  »*2*CC( 1 , 1 ) 

1 +2. * ( AA ( 2, 1 ) -AA (1,1 ) )*SUM1+SUM2+2.»SUM3) 

DRAFT  =  YCSO( I  I ) -YCL < 1  I ) 

AHCYL=DENS I TY»P I *DRAFT**2/2 . 

CH=AH/AHCYL 
PRINT  288, DRAFT 

288  FORMAT ( *0*9X#SECT I  ON  DRAFT  =*F7.2) 

GO  TO  292 

290  AV1=DENSITY*(2.*PI*AA(1,1  )  *  <  AA  ( 1  , 1  )  +AA  (3,1  ))-AREA(  ID) 
AV2=DENSITY»(2.*PI*AA(1 ,1 )*(AA(1,1 )+AA( 3, D ) -AREAMF ) 

AVCYL=DENS I TY*P I *BD2**2 

CV=AV1 /AVCYL 

AH=DENS  I  TY*(2.  *PI  *AA(  1 , 1  )  *  ( AA  ( 1  ,1  )-AA(3, 1  )  )-AREA(  ID) 

H=DCU( I  I ) -YCL ( I  I  ) 

HD2=H/2 . 

AHCYL=DENS I TY*PI *HD2**2 
CH=AH/ AHCYL 
PRINT  291, H 

291  FORMAT ( *0*9X#SECT I  ON  HEIGHT  =»F7.2) 

292  PRINT  293, AVI, AH 

293  FORMAT ( »0*T44»VERT I CAL*T59*HOR I ZONTAL*/ 
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NDATA( 5) =0 

Q 

C***  CALCULATE  CX  AND  CY  FUNCTIONS  OF  ANGLES  IN  ZETA  PLANE 
DO  220  M=1 ,MM 

THETAM=THSTP*FLOAT (M-1 )-PID2 
DO  220  N=1 ,NN 

ARG= ( ALPHA -BETA*FLOAT (N) ) * THE TAM 
CARG=COS ( ARG ) 

SARG=SIN(ARG) 

I F ( SUBMERG )  GO  TO  214 
CX(N,M)=CARG 
CY(N,M)=SARG 
GO  TO  220 

214  I F ( MOD (N,2).EQ.O)  GO  TO  216 
CX(N,M)=CARG 
CY(N,M)=SARG 
GO  TO  220 
216  CX(N,M)=-SARG 
CY(N,M)=CARG 
220  CONTINUE 

Q 

C***  CALCULATE  MAPPING  FUNCTION  PARAMETERS  USING  LS  CRITERION 
N=NN 

IF(NN.LE.NNDFALT)  GO  TO  224 
N=NNDFALT 
IF(IN.NE.O)  N= I N 
224  NDATA(6)=NDA( I  I ) 

230  NDATAC 1 )=N 
240  DO  242  K=1,NN 
AA ( K , 1 )=AI (K) 

242  CONTINUE 

CALL  LSQSUB(NDATA,X,AA,LSQFUN,1 .E-5) 

IF ( AA ( 26, 1 ) .NE.O. )  GO  TO  244 
IF(NDATA(6).EQ.1 )  GO  TO  250 
NDATA(6)=1 
GO  TO  240 

244  PRINT  245,AA(26,1),NDATA(6),(AA(I,1),AA(I,2),I=1,N) 

245  FORMAT ( 1  HO,  9X*FINAL  MAPPING  PARAMETERS  AFTER  *F3.0*  ITERATIONS 
+(  QME  )*10X*M0DE*I2/(7X,5(3X,F11.7*(*F8.6*)*))) 

I F (A  I  PR  I  NT)  PR  I  NT (3, 1 05)  N ,  ( AA  ( I  ,  1  ) ,  I  =  1 ,  N ) 

N=N+1 

IF(N.LE.NN)  230,256 

250  PRINT  255,NDATA(4 ) , (AA( I , 1 ) , AA ( I , 2 ) , I =1 ,NN) 

255  FORMAT (1  HO, 1X*N0  CONVERGENCE.  MAPPING  PARAMETERS  AFTER* 1 3*  ITERATI 
+ONS  (  QME  )*/ ( 7X, 5( 3X,F 1 1 . 7*(*F8.6*)* ) ) ) 

256  PRINT  257,AA(26,2) 

257  FORMAT (1  HO, 9X*QME  OF  OVERALL  FIT*F10.7) 

C#*»  PRINT  CALCULATED  SHIP  FORM  IN  COORDINATES  RELATIVE  TO  SHIP  ORIGIN  (SO) 
PRINT  260 

260  FORMAT (*Q*9X*CALCULATED  SHIP  F0RM*/14X,4(*X*6X*Y*5X*RM-RS*7X) ) 

DO  262  M=1 ,MM 
X(8,M)=X(9,M)-X(2,M) 

262  CONTINUE 

PRINT  264, (X(6,M),X(7,M),X(8,M),M=1 ,MM) 
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110  NN=NNA( I  I ) 

GO  TO  120 
115  NN=NNDFALT 
120  DO  130  J=1,JJ 
P=PC I  N(  J, I  I ) 

CALL  SPLCAL 
D I F  ( J ) =RC I N ( J ,  I  I  )-R 
130  CONTINUE 

PRINT  135,JJ,KK,(DIF(J),J=1,JJ) 

135  FORMAT ( 1  HO, 9X*RES I  DUALS  FOR  SPLINE  FIT  TO* I  3*  INPUT  POINTS  WITH* I  2 
+  *  POINTS  PER  ARC  OF  SPL I NE*/ ( 1  OX, 1 0G1 2. 3 ) ) 

PRINT  1 40, E (II) 

140  FORMAT (1  HO, 9X*RMS  ERROR  OF  SPLINE  FIT  =*G12.4) 

C 

C*»*  PREPARE  COEFFICIENTS  AND  CONSTANTS  FOR  FITTING  MAPPING  FUNCTION 
MM=2 . 5*F  LOAT ( NN ) 

IF(MM.LT.JJ)  MM= J J 
THSTP=PI /FLOAT (MM-1 ) 

I F ( SUBMERG)  GO  TO  160 
THSTP=THSTP/2. 

150  ALPHA=3. 

BETA=2. 

IF(IN.NE.O)  GO  TO  170 
XL IM=DCU( I  I ) 

BD2=XL IM-XCSO( I  I ) 

D=YCSO( I  I )-YCL(  I  I  ) 

A I (1 )=(D+BD2)/2. 

A I (2)=(D-BD2)/2. 

JN=3 

GO  TO  170 
160  ALPHA=2. 

BETA=1 . 

CALL  SUBBEAM (BD2) 

IF( IN.NE.O)  GO  TO  170 
YSEMI  =  (DCU( I  I ) - YCL ( I  I ))/2. 

A I (1 )  =  (YSEMI+AREA( I  I ) /PI /YSEMI )/2. 

A I ( 2 )  =  ( DCU ( I  I )+YCL(l I ) )/2. -YCSO( I  I ) 

A I ( 3 ) =A I ( 1 ) -YSEMI 
JN=4 

170  IF ( IN.NE.O)  JN= I N+1 
IF ( JN.GT.NN)  GO  TO  185 
DO  180  N  =  JN,NN 
A I (N ) =0. 

180  CONTINUE 
185  DO  190  N=1 ,NN 
AA(N,3)=N 
190  CONTINUE 

PRINT  200 , ( A I (N) ,N=1 , NN ) 

200  FORMAT (1  HO, 9X* INITIAL  MAPPING  PARAMETERS*/ ( 7X, 5( 3X,F 1 1 . 7, 1 0X) ) ) 

DO  210  M= 1 ,MM 
X( 1 ,M)=M 
210  CONTINUE 
NDATA( 2 ) =MM 
NDATA(3)=2 
NDATA(4)=50 
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JJ=J JA( I ) 

E( I )=0. 

XCSO( I )=SGN( I )*XCS0( I ) 

KKA( I )={ JJA( I )+22)/23+1 
I F (KKA( I ) .LT.3)  KKA( I )=3 
DO  40  J=1,JJ 

IF (SGN( I ) .EQ.-1 )  PC  I N ( J , I ) =  S I GN (PI, PC  I N ( J , I ) )-PC I N ( J , I ) 

IF(J.EQ.I)  GO  TO  40 

IF(PCIN(J, I ) .GT. PC  I N( J-1 , I ))  GO  TO  40 

PRINT  35,  NS( I ) , J 

35  FORMAT (*0+++++++  IMPROPER  DATA  FORMAT  OR  CALCULATION  ORIGIN  FOR  SE 
+CTION*l 5*  NEAR  POINT  NUMBER*! 5*  ...  STOP  ++++++++++*) 

STOPPER  =1 . 

40  CONTINUE 
50  CONTINUE 

IF(STOPPER.EQ. 1 . )  STOP  °BY  SHIPFIT0 
IF (PLOT)  CALL  SH I  PLOT ( 1 ) 

C 

C***  FIT  CUBIC  SPLINE  FUNCTIONS  TO  CO  REFERENCED  DATA 
DO  70  1=1,111 
J J= J JA( I ) 

CALL  SPLSQ1 ( J  J , PC  I N ( 1 , I ) , RC I N ( 1 , I ) ,CF ( 1 , I ) ,E ( I ) ,KKA( I ) ) 

DO  60  J=1,JJ 
RC2( J, I ) =RC I N ( J , I )**2 
60  CONTINUE 

CALL  SPLSQ1 ( J  J , PC  I N ( 1 , I ) ,RC2( 1 , I ) ,CA( 1 , I ) ,-1 . ,KKA( I ) ) 

70  CONTINUE 
C 

C***  CALCULATE  SUBMERGED  CROSS-SECTIONAL  AREAS  AND  SET  LOCATION  OF  SHIP  ORIGIN. 
C  IF  VESSEL  IS  FLOATING  FIND  EQUILIBRIUM  POSITION.  OBTAIN  LIMITS. 

CALL  ORIENT 
C 

C***  CALCULATE  ADDED  MASSES  OF  SECTIONS 

READ  75, NCASES , ( NSP I CK ( J ) , J  =  1 , NCASES ) 

75  F0RMATU6I5) 

DO  300  IS=1, I  I  I 
I  I  =  I S 

IF (NCASES. EQ.O)  GO  TO  90 
DO  80  J=1, NCASES 
I F ( NS (II) ,EQ.NSPICK( J ) )  GO  TO  90 
80  CONTINUE 
GO  TO  300 
90  JJ=JJA( I  I ) 

KK=KKA( I  I ) 

XSHIFT=XCSO( I  I  ) 

YSH I FT=YCSO( I  I  ) 

PRINT  95, (Tl TLE( J) , J=1 ,8) ,NS ( I  I ), FRAME ( I  I ) 

95  FORMAT (*1*4X8A10//5X*CALCULATION  OF  ADDED  MASSES*// 

1  5X*STATI0N  NUMBER* 15,1 0X*L0CATED  AT  FRAME*F7.2) 

I  N=0 

IF (NNA( I  I ) ) 100, 115,110 
100  NN=-NNA( I  I ) 

READ  105, IN, (Al(l), 1=1, IN) 

105  FORMAT ( 1 1 0/ ( 8G 10.4) ) 

GO  TO  120 
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PROGRAM  SH I PF ! T( INPUT, OUTPUT, DATA1 ,DATA2,AIF ILE, 

1  TAPE1 =DATA1 ,TAPE2=DATA2,TAPE3=AIF ILE) 

DIMENSION  NDAT A ( 6 ) , NNA ( 24 ) , SGN (24),DIF(64),NSPICK(24), 

1  A I (24) ,NDA(24) 

COMMON  /BL0CK1  /  SUBMERGES  (24 ),  FRAME  ( 24 )  ,KKA(  24 )  ,KK,  JJA(24 ),  JJ , 

1  XCS0(24) , YCS0(24) ,RC IN(64,24) ,PC IN<64,24) ,CF(28,24) ,E(24) ,DCU(24) 
2, YCL(24),AREA(24), DENSITY, P,R,RP,RPP, 11,111 
COMMON  /BL0CK2/  AA( 26, 3) ,TI TLE( 8 ) , XSH I  FT, YSH I  FT, ALPHA, BETA, 

1  CX(24,64 ) ,CY(24, 64 ) ,MM, NN.NSEGS 

COMMON  /BL0CK3/  PID1 80, YDC0(24) ,CA ( 28, 24 ) ,RC2(64,24) ,X(9, 64) 

EXTERNAL  LSQFUN 

LOGICAL  SUBMERG, A I  PR  I  NT , PLOT 

DATA  PI ,PID2/3. 141 5926535898,1 .5707963267949/ 

P I D 1 80=P 1/180. 

READ  5, (T I TLE ( I ) , I =1 , 8 ) 

5  FORMAT (8A1 0) 

READ  6, SUBMERG, NNDFALT, PLOT, A I  PR  I  NT, SPVOL 

6  FORMAT(L10, I  1 0, 2L1 0,E1 0. 0) 

C 

C*  SET  DEFAULT  VALUE  OF  NNDFALT. 

IF(NNDFALT.EQ.O)  NNDFALT=8 

C*  SET  DEFAULT  VALUE  OF  SPVOL  TO  35  CU  FT  /  TON  (LONG) 

IF ( SPVOL. EQ.O. )  SPVOL=35. 

DENS  I TY=1 ./SPVOL 
C 

C***  READ  AND  CONVERT  X,Y  SHIP  FORM  DATA  FROM  DATA  ORIGIN  (DO) 

C  TO  CALCULATION  ORIGIN  (CO) 

111=0 

10  READ ( 1 ,*)  I D I G , XD I G , YD  I G 
I F  ( EOF  (1  ) .  NE .  0 )  GO  TO  30 
IF( IDIG.NE.1 )  GO  TO  20 
111=111+1 

READ  1 5,NS( I  I  I ), FRAME ( I  I  I ) , SCALE, NNA ( I  I  I ) ,NDA( III) 

15  FORMAT (110, 2E 10. 0,2110) 

C*  SET  DEFAULT  VALUE  OF  SCALE  TO  1. 

IF (SCALE. EQ.O. )  SCALE=1 . 

JJA( I  I  I )=0 
XDM I N=1 . E 1 0 
XDMAX=-1 . E 1 0 
XDCO=XD I G*SCALE 
YDCO( I  I  I )=YD IG*SCAL£ 

XCSO( I  I  I )=-XDCO 
GO  TO  10 

20  J J=J JA( I  I  I )=J JA( I  I  I )  +  1 
XD=XDIG*SCALE 
XDMIN=AMIN1 ( XD , XDM I N ) 

XDMAX=AMAX1 (XD,XDMAX) 

SGN( I  I  I ) =S I GN ( 1 . , XDM 1 N+XDMAX ) 

XC=XD-XDCO 

YC=YD I G*SCALE -YDCO ( III) 

RC I N ( J  J , I  I  I )=SQRT(XC**2+YC*»2) 

PCIN( JJ, I  I  I )=ATAN2( YC,XC) 

GO  TO  10 
30  ST0PPER=0. 

DO  50  1=1,111 
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however,  there  is  no  special  point  allocation  feature  of  this  kind  built  into 
SHIPFIT.  Satisfactory  results  can  usually  be  achieved  by  simply  increasing  the 
number  of  mapping  function  coefficients  desired,  since  the  number  of  image 
points  is  increased  proportionately. 

For  all  difficult  sections,  plotting  of  the  input  points  and  mapping 
function  is  recommended.  Another  useful  check  of  the  level  of  success  achieved 
in  fitting  the  mapping  function  to  the  spline  function  is  provided  by  a  "check" 
value  of  the  vertical  added  mass.  This  is  printed  automatically  for  all  cross 
sections.  The  "check"  value  is  obtained  by  using  a  cross-sectional  area 
calculated  from  the  full  set  of  fitted  mapping  coefficients,  while  the  value 
labeled  "best"  is  obtained  by  using  an  area  calculated  from  the  cubic  spline 
function.  A  discrepancy  between  the  "best"  and  "check"  values  indicates  a 
problem  in  fitting  the  mapping  function  coefficients.  Most  sections  will  be 
fitted  easily,  without  the  need  for  the  special  features  built  into  SHIPFIT. 
Additional  form  data  that  may  be  used  as  test  cases  for  the  program  can  be  found 
in  Reference  1. 


FIGURE  11.  PLOT  OF  INPUT  POINTS  AND  SPLINE  FUNCTION  PRODUCED  BY  INPUT  DATA  OF  FIGURE  7 
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FIGURE  10.  CALCULATED  SHIP  ORIENTATION  FOR  SAMPLE  PROBLEM  DATA 
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SUBROUTINE  AREACAL 

C***  AREACAL  CALCULATES  THE  SUBMERGED  CROSS-SECTIONAL  AREA  USING  THE  POLAR 
C  COORDINATE  SYSTEM  CENTERED  AT  THE  CALCULATION  ORIGIN  (CO). 

COMMON  /BLOCK 1 /  SUBMERG, NS(24),FRAME(24),KKA(24),KK, JJA(24), JJ, 

1  XCSO ( 24 ) , YCSO ( 24 ) , RC I N ( 64 , 24 ) , PC  I N ( 64 , 24 ) , CF ( 28 , 24 ) , E ( 24 ) , DCU ( 24 ) 

2, YCL(24),AREA(24), DENSITY, P,R,RP,RPP, 11,111 
COMMON  /BLOCK3/  PID180,YDCO(24) ,CA(28,24) ,RC2(64,24) ,X(9,64) 

LOGICAL  SUBMERG 

DATA  PI ,PID2/3. 1415926535898, 1.5707963267949/ 

IF (SUBMERG)  10,30 
10  AY=0. 

I F ( XCSO ( I  I J.NE.O. )  GO  TO  20 

P=PCL=-PID2 

CALL  SPLCAL 

YCL( I  I )=-R 

P=PCU=PI D2 

CALL  SPLCAL 

DCU( I  I )=R 

AX=0. 

GO  TO  80 

20  CALL  YL I M I T ( PC  I N ( 1 , I  I ) , YCL (II)) 

CALL  YL I M I T ( PC  I N ( J  J , I  I ) , DCU (II)) 

PCL=ATAN2  ( YCL  (  I  I  ) ,  XCSO  (ID) 

PCU= ATAN2  ( DCU  (  I  I  ) ,  XCSO  (ID) 

AX=XCSO( I  I )*(DCU( I  I ) —YCL ( I  I) )/2. 

GO  TO  80 

30  I F ( XCSO ( I  D.NE.O.)  GO  TO  40 
P=PCL=-PID2 
CALL  SPLCAL 
YCL( I  I ) =-R 
AX=0. 

GO  TO  50 

40  CALL  YL  I M I T  ( PC  I N  ( 1 ,  I  I) ,  YCL  (ID) 

PCL=ATAN2(YCL(  I  I  )  ,XCSO(  ID) 

AX=XCSO( I  I )*( YCSO( I  I )-YCL( I  I ) )/2. 

50  IF(YCSO(  I  D.NE.O.)  GO  TO  60 
P=PCU=0. 

CALL  SPLCAL 
DCU( I  I )=R 
AY=0. 

GO  TO  80 
60  DO  65  J=1,JJ 

I F ( RC 1 N ( J , I l)*SIN(PCIN(J, I l))-YCSO( I D.GT.O. )  GO  TO  70 
65  CONTINUE 
J  =  J  J 

70  CALL  XL  I M I T  ( PC  I N  ( J ,  I  I  ) ,  DCU  (ID) 

PCU=ATAN2 ( YCSO ( I  I) , DCU ( I  I) ) 

AY=YCSO( I  I )*(XCSO( I  I ) -DCU ( I  I) )/2. 

80  CALL  SPL I  NT ( J  J , PC  I N ( 1 , I  I ) ,RC2 (1,1  I ) ,KKA( I  I ) , AC02,PCL,PCU,CA( 1 , I  I ) ) 
AREA( I  I ) =2 . * ( AC02/2 . -AX-AY ) 

RETURN 

END 
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FUNCTION  CC( IU, IV) 

IR=2* IU-1 
I S=2* I V-1 

IFdR.EQ.  IS)  GO  TO  30 
IF(IR.LT.IS)  GO  TO  10 
I  X=  I  s 
I  S=  I R 
I  R=  I X 

10  FNF=NF= I S- IR 
SUM=0. 

DO  20  N=2,NF,2 

SUM=SUM+1 ./FLOAT ( IR+N)/FLOAT( IS-N) 
20  CONTINUE 

CC=-IR*IS/2./FNF*SUM  ' 

GO  TO  50 
30  NF=2* I S 
SUM=0. 

DO  40  N=2,NF, 2 
SUM=SUM+1 ./(FLOAT! IS-N) )**2 
40  CONTINUE 
CC= I  S/4. *SUM 
50  RETURN 
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SUBROUTINE  LSQFUN (NDATA, X, D1 ,A,D2) 

C***  LSQFUN  CALCULATES  THE  MAPPING  FUNCTION  AND  SPLINE  FUNCTION  RADII  (RELATIVE 
C  TO  THE  CALCULATION  ORIGIN)  AND  THE  DERIVATIVES  OF  THE  RESIDUAL  TERM  AS 
C  REQUIRED  BY  LSQSUB. 

COMMON  /BLOCK  1  /  SUBMERGES  (24 )  .FRAME  (24 )  ,KKA(24 )  ,KK,  JJA(24 ),  JJ, 

1  XCSO ( 24 ) , YCSO ( 24 ) , RC I N ( 64 , 24 ) , PC  I N ( 64 , 24 ) , CF ( 28 , 24 ) , E ( 24  > , DCU ( 24 ) 

2, YCL(24),AREA(24), DENSITY, P,R,RP,RPP, 11,111 
COMMON  /BL0CK2/  AA(26, 3) , Tl TLE( 8 ) , XSH I  FT, YSHI FT, ALPHA, BETA, 

1  CX(24,64 ) ,CY(24,64 ) ,NN 

DIMENSION  NDATA(6),X(9),D1 (1 ),A( I ) ,DR(24 ) ,DP(24 ) ,D2(28,28) 

N1 =NDATA( 1 ) 

M=X(  1  ) 

XM=YM=0. 

DO  10  N=1 ,N1 
XM=XM+A ( N ) *CX ( N , M ) 

YM= YM+A ( N ) *CY ( N ,M ) 

1 0  CONT I NUE 
X ( 6 ) =XM 
X(7)=YM 
XXC=XM+XSH I  FT 
YYC=YM+YSHI FT 
RC2=XXC**2+YYC**2 
RC=SQRT(RC2) 

CXI =XXC/RC 
CY1 =YYC/RC 
CX2=CX1 /RC 
CY2=CY1/RC 
P=ATAN2 ( YYC, XXC ) 

CALL  SPLCAL 
X(2)=R 
X ( 9 ) =RC 
DO  20  N=  1 ,  N 1 

DR(N)=CX1 #CX(N,M)+CY1#CY(N,M) 

DP ( N ) =CX2*CY ( N , M ) -CY2*CX ( N , M ) 

D1 ( N ) =DR ( N ) -RP*DP ( N ) 

20  CONTINUE 

I F (NDATA(6) )25,40 
25  DO  30  L= 1 , N 1 
DO  30  K=L,N1 

D2R=(CX(L,M)*CX(K,M)+CY(L,M)*CY(K,M)-DR(K)*DR(L))/RC 

D2P= ( CX ( L , M ) *CY ( K , M ) -CY ( L , M ) *CX ( K , M ) ) /RC2-2 . *DR ( L ) *DP ( K ) /RC 

D2(K,L) =D2R-RP*D2P-RPP*DP (K)*DP(L) 

30  CONTINUE 
40  RETURN 
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SUBROUTINE  LSQSUB  (N,  X,  A,  VF,  DELT) 

0 

C***  LSQSUB  IS  A  GENERAL  NON-LINEAR  LEAST  SQUARES  FITTING  ROUTINE. 

C  IT  WAS  WRITTEN  BY  TED  ORLOW  AND  IS  AVAILABLE  THROUGH  THE  NSWC/WHITE 
C  OAK  PROGRAM  LIBRARY. 

C 

DIMENSION  DF (25 ) ,  A(26,3),  X(9,1),  V(28,28),  N(5),  D2F(28,28) 

NPAR  =  N(1 ) 

NPARM1 =NPAR-1 
NPT  =  N (2) 

NPCD  =  N ( 3 ) 

NWT  =  NPCD+1 
NWAR  =  N(5) 

NIT  =  0 

I F ( NP AR . LT . 26 . AND . NPAR . GT . 0 . AND . NPCD . LT . 8 . AND . NPCD . GT . 1 ) GO  TO  20 
PRINT  10,  NPAR, NPCD 

10  FORMAT (1  HO, 6HNPAR  =, I  5, 4X, 2H0R, 2X, 6HNPCD  =,I5,4X, 

115HIS  OUT  OF  RANGE) 

STOP 

20  IF (NWAR)  70,30,50 
C***  ABSOLUTE  FIT 

30  DO  40  K  =  1,  NPT 
40  X(NWT,  K)  =  1. 

SUMWT=NPT 

C=1. 

GO  TO  100 
C***  RELATIVE  FIT 

50  DO  60  K  =  1,  NPT 
60  X(NWT,  K  )  =  1./  X(NPCD,  K  )  **2 
C***  WEIGHTS  SUMMED 
70  SUMWT=X(NWT , 1 ) 

DO  80  L  =  2,  NPT 
80  SUMWT  =  SUMWT  +  X(NWT,  L) 

C=SQRT ( SUMWT/FLOAT ( NPT ) ) 

C**»  CALCULATE  QME  OF  FIT  BEFORE  (A(26,3))  AND  AFTER  (A(26,2>>  ITERATING 
100  SUM2  =  0. 

DO  110  1=1, NPT 

CALL  VF  (  N  ,  X(1,l),  DF,  A  ,D2F) 

110  SUM2  =  SUM2  +  X(NWT, I )  *  (X(9, I >-X(NPCD, I ))**2 
A ( 26 , 2 ) = SQRT ( SUM2/ SUMWT ) 

I F (N(4 ) . EQ.O)  GO  TO  220 
IF(NIT.NE.O)  GO  TO  190 
NIT  =  N(4) 

A (26, 3)  =  A ( 26 , 2 ) 

C***  MAJOR  ITER  LOOP 

DO  180  K  =  1,  NIT 

C**»  SET  UP  SYSTEM  OF  LINEAR  EQUATIONS. . A( I ,2 )-V( I , J )*( A( J, 1 )-A°( J, 1 ) )=0 
DO  120  J  -  1,  NPAR 
A( J , 2)  =  0. 

DO  120  I  =  J , NPAR 
D2F ( I , J )=0. 

120  V(l,J)=0. 

DO  150  KK=1,NPT 

CALL  VF  (  N  ,  X ( 1 , KK ) ,  DF,  A  ,D2F) 
WDELY=(X(9,KK)-X(NPCD,KK))*X(NWT,KK) 


ooo  o  o  ooo  ooo 
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DIVIDE  PIVOT  ROW  BY  PIVOT  ELEMENT 

310  PIVOT  =A( ICOLUM, ICOLUM) 
DETERM=DETERM#P I VOT 
330  A( ICOLUM, ICOLUM)=1 .0 
340  DO  352  L=1,N 

IF  (PIVOT  .NE.  0.  )  GO  TO  350 
A( I COLUM,  L)  =  0. 

GO  TO  352 

350  A(l COLUM, L)=A( I COLUM, L)/PI VOT 

352  CONTINUE 

355  I F (M)  380,  380,  360 

360  DO  375  L=1,M 

IF  (PIVOT  .NE.  0.  )  GO  TO  370 
B( ICOLUM, L)  =  0. 

GO  TO  375 

370  B( ICOLUM, L)=B( ICOLUM, L ) /PI VOT 
375  CONTINUE 

REDUCE  NON-PIVOT  ROWS 

BMAX=0. 0 

380  DO  550  LI =1 ,N 

390  IF  (LI  .EQ.  ICOLUM)  GO  TO  550 

400  T=A( LI, ICOLUM) 

420  A(L1 , I COLUM) =0.0 
430  DO  450  L=1  ,N 

SUB  =A( ICOLUM, L)*T 
A(L1,L)=A(L1,L) -SUB 


IF( I NDEX (LI ,3) .EQ. 1 )  GO  TO  450 
IF(ABS(A(L1 ,L))-EPS4*  ABS(SUB) ) 449, 449, 450 

449  BMAX=AMAX 1 ( BMAX , ABS ( A ( L 1 , L ) ) ) 

450  CONTINUE 
455  I F (M)  550,  550,  460 
460  DO  500  L=1,M 

500  B(L1,L)=B(L1,L)-B( ICOLUM, L)*T 
550  CONTINUE 

INTERCHANGE  COLUMNS 

600  DO  710  1=1, N 
610  L=N+1 - 1 

620  IF  ( I NDEX( L , 1 )  .EQ.  I NDEX(L, 2 ) )  GO  TO  710 
630  JR0W= I NDEX(L, 1 ) 

640  JC0LUM= I NDEX ( L , 2 ) 

650  DO  705  K=1,N 
660  SWAP=A(K, JROW) 

\  670  A ( K , JROW ) =A ( K , JCOLUM ) 

:  700  A ( K, JCOLUM )= SWAP 

705  CONT I NUE 
I-  710  CONTINUE 

•  DO  730  K  =  1 ,N 
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IF  <INDEX(K,3)  .NE.  1)  GO  TO  715 
730  CONTINUE 
RETURN 
715  1 0=2 
RETURN 
END 


A- 13 
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SUBROUTINE  MATRI Xf A,N1 , 11 ,8, Ml , DETERM, ID) 

C 

C»**  MATRIX,  WRITTEN  BY  CHARLES  NEWMAN,  IS  AVAILABLE  ON  THE  NSWC/WHITE  OAK 
PROGRAM  LIBRARY. 

MATRIX  INVERSION  WITH  ACCOMPANYING  SOLUTION  OF  LINEAR  EQUATIONS 

TEST  FOR  LOST  OF  DIGITS  DUE  TO  SUBTRACTION 

DIMENSION  A( I  1 , 1 1 ),B( 11 ,M1 ) , I NDEX( 50, 3 ) 

EQUIVALENCE  (IROW,JROW),  ( ICOLUM, JCOLUM) ,  (AMAX,  T,  SWAP) 

DATA  EPS4/1.E-7/ 

INITIALIZATION 


BMAX=0.0 

10=1 

M=M1 

N=N1 

10  DETERM* 1.0 
15  DO  20  J= 1 , N 
20  INDEX( J , 3 )  =  0 
30  DO  550  1=1 ,N 

SEARCH  FOR  PIVOT  ELEMENT 

40  AMAX=0. 0 
45  DO  105  J  —  1 , N 

IF  ( I NDEX( J , 3 )  .EQ.  1 )  GO  TO  105 
60  DO  100  K=1 ,N 

IF ( INDEX(K,3)-1 )  80,  100,  100 
80  IF  (AMAX  ,GE.  ABS(A(J,K)))  GO  TO  100 
85  I ROW=  J 
90  IC0LUM=K 

AMAX  =  ABS ( A ( J , K ) ) 

100  CONTINUE 
105  CONTINUE 

INDEX( ICOLUM, 3)  =  INDEX( ICOLUM, 3)  +1 
260  INDEX( I , 1 )= IROW 
270  I NDEXC I ,2)= ICOLUM 
IF (BMAX.EQ.AMAX) ID=2 

INTERCHANGE  ROWS  TO  PUT  PIVOT  ELEMENT  ON  DIAGONAL 

130  IF  (IROW  .EQ.  ICOLUM)  GO  TO  310 
140  DETERM=-DETERM 
150  DO  200  L=1 ,N 
160  SWAP=A( IR0W,L) 

170  A( IROW,L)=A( ICOLUM, L) 

200  A( ICOLUM, L)=SWAP 
I F (M)  310,  310,  210 
210  DO  250  L=1 ,  M 
220  SWAP=B( IROW,L) 

230  B ( I ROW, L ) =B ( I COLUM, L ) 

250  B( ICOLUM, L )=SWAP 
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C***  COMPUTE  A( I ,2)  VECTOR  AND  V  MATRIX  ELEMENTS 
DO  130  J=1 ,  NPAR 
A(J,2)=A(J,2)-WDELY*DF(J) 

DO  130  I  =  J , NPAR 

V( I ,J)=V( l,J)+X(NWT,KK)*DF( I ) *DF ( J ) +WDEL Y*D2F ( I ,  J  ) 

130  CONTINUE 
150  CONTINUE 

DO  160  J=1 , NPARM1 
DO  160  I =J,NPARM1 
160  V(J, 1+1 )=V(I+1,J) 

C***  INVERT  MATRIX  AND  COMPUTE  CHANGE  IN  °A°  VECTOR 

CALL  MATRIX  {  V,  NPAR,  28,  A{1, 2),  1,  DETERM,  ID) 

C***  CALCULATE  IMPROVED  °A°  VECTOR  AND  TEST  FOR  CONVERGENCE 
C***  IF  NO  CONVERGENCE  A(26,1)=0  ...  IF  CONVERGENCE  AC 26, 1 ) =K 
A  ( 26 , 1 )  =  K 
DO  170  I  =  1,  NPAR 
M  =  A(l, 3) 

A(M, 1 )  =  A(M, 1 )  +  A( I ,2) 

IF  (A(26, 1 )  .E<?.  0.)  GO  TO  170 
SSS=ABS ( A( I ,2)) 

IF(A(M, 1 ) .NE.O. )  SSS=AMI N1 (A8S(A( I ,2)/A(M, 1 ) ) ,SSS) 
IF(SSS.GE.DELT)  A(26,1 )=0. 

170  CONTINUE 

IF  (AC 26 , 1 )  .NE.  0.)  GO  TO  100 
180  CONTINUE 
GO  TO  100 

C**»  Ad, 2)  NOW  BECOMES  QME  OF  I -TH  PARAMETER 
190  DO  210  1=1, NPAR 
M  =  A(  I  ,  3  ) 

I F( V ( I , I )  .GT.  0.)  GO  TO  200 
A(M,2>  — 1. 

GO  TO  210 

200  ACM, 2)  =  A ( 26 , 2 )  *  SQRT ( V ( I , f  7  > 

210  CONTINUE 
220  RETURN 
END 
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SUBROUTINE  ORIENT 
DIMENSION  WA ( 3 , 48 ) , I WA ( 3 , 48 ) 

EQUIVALENCE  (WA.IWA) 

COMMON  /BLOCK  1  /  SUBMERGES  (24 )  .FRAME  ( 24 )  ,KKA  (24 )  ,KK,  JJA  (24 ),  JJ, 

1  XCS0(24) ,YCS0(24) ,RC IN(64,24) ,PC IN(64,24) ,CF (28,24) ,E(24) ,DCU(24) 

2, YCL(24),AREA(24), DENSITY, P,R,RP,RPP, 11,111 
COMMON  /BL0CK2/  AA ( 26,3) ,TI TLE(8 ) ,XSHIFT,YSHI FT, ALPHA, BETA, 

1  CX(24,64 ) ,CY(24,64 ) ,MM,NN,NSEGS 

COMMON  /BL0CK3/  P ID1 80, YDCO( 24) ,CA( 28, 24 ) , RC2(64, 24 ) ,X(9,64 ) 

LOGICAL  SUBMERG,EQCAL, FLOODED 
REAL  LBP 

IF (SUBMERG)  10,30 
C 

C***  VESSEL  IS  TOTALLY  SUBMERGED. 

C  CALCULATE  TOTAL  SECTION  AREAS  AND  LOCATE  THE  SHIP  ORIGIN  AT  THE 
C  GEOMETRIC  CENTERS. 

10  DO  20  I S=1 , I  I  I 

I  I  =  I S 

J J=JJA( I  I ) 

KK=KKA( I  I  ) 

CALL  AREACAL 

YCSO( I  I )  =  ( DCU ( I  I )+YCL( I  I ))/2. 

20  CONTINUE 
RETURN 
C 

C**»  VESSEL  IS  FLOATING  AT  SURFACE. 

30  READ  35, LBP, NSEGS, EQCAL 
35  FORMAT (E10.0, 110, L10) 

C*  SET  DEFAULT  VALUE  FOR  NSEGS 
I F (NSEGS. EQ.O)  NSEGS=20 
SEGW I D=LBP/FLOAT ( NSEGS ) 

SEGDEN=SEGW I D*DENS I TY 
IF (EQCAL)  GO  TO  100 
C 

C*»*  SHIP  EQUILIBRIUM  POSITION  IS  NOT  CALCULATED...  INPUT  WATERLINE  IS  USED. 

C  CALCULATE  ASSOCIATED  BUOYANCY,  FREEBOARD,  DRAFT  FOR  EACH  SECTION  AND  PRINT 
READ  40,WATERLN 
40  FORMAT(EIO.O) 

BSUM=0. 

DO  50  1=1,111 

YCSO( I )=WATERLN-YDCO( I ) 

J  J = J  J  A  (  I  ) 

KK=KKA ( I ) 

II  =  1 

CALL  AREACAL 

C*  WA  AND  I WA  ARRAYS  ARE  USED  FOR  SORTING/PRINTING  RESULTS. 

WA ( 1 , I ) =F  LOAT ( NS ( I ))»SEGWID 
WA(2, I )=AREA( I )*$EGDEN 
I WA( 3, I )=l 
BSUM=BSUM+WA ( 2 ,  I ) 

50  CONTINUE 

CALL  SORT (WA, I  I  1*3,3) 

PRINT  60, ( T I TLE ( J ) , J= 1 ,8) 

60  F ORMAT ( *  1 *4X , 8A 1 0// 5X* SH I P  EQUILIBRIUM  POSITION  NOT  CALCULATED  ... 

1  PROGRAM  USES  INPUT  WATERLINE*) 
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PRINT  65,WATERLN,BSUM 

65  FORMAT («0*9X* INPUT  WATERLINE  (RELATIVE  TO  DATA  ORIGIN)  =*F6.2 
1 /*0*9X*T0TAL  BUOYANCY  FORCE  =*F10.2) 

PRINT  70 

70  FORMAT ( *0*// 1 4X*STAT I 0N*3X*FRAME*3X*FREEB0ARD*3X*DRAFT*3X 
1*0 1  ST  FROM  STA  0*3X*BU0YANCY*) 

DO  90  1=1,111 
K= I WA( 3, I ) 

J=JJA(K) 

DRAFT=YCSO ( K ) -YCL ( K ) 

FREE8D=RC IN(J,K)*SIN(PCIN(J,K)) -YCSO(K) 

PRINT  80,NS(K) , FRAME (K ) ,FREEBD, DRAFT, WA( 1 , I ) ,WA(2, I ) 

80  FORMAT ( 16X I 2,6XF5. 1 , 5XF5. 2, 5XF5. 2, 7XF6. 2,8XF8. 2 ) 

90  CONTINUE 
GO  TO  300 
C 

C**»  FIND  EQUILIBRIUM  POSITION  OF  SHIP  BY  BALANCING  SUMS  OF  WEIGHT  AND 
C  BUOYANCY  FORCES  (WSUM  AND  BSUM)  AND  MOMENTS  (WMOM  AND  BMOM). 

C  CALCULATE  ASSOCIATED  BUOYANCY,  FREEBOARD,  DRAFT  FOR  EACH  SECTION  AND  PRINT 
100  LL=0 

READ  110,  COLWID,WSCALE,ZSCALE, FLOODED 
110  FORMAT(3E10.0,L10) 

I F (COLW ID.GT. 1.01 *SEGW ID)  PRINT  115 

115  FORMAT (*...  WARNING  WEIGHT  HISTOGRAM  COLUMN  WIDTHS  EXCEED  THE  SHIP 
1  SEGMENT  LENGTHS.  ORIENT  ASSUMES  THE  CONTRARY  ...*) 

C*  SET  DEFAULT  VALUE  FOR  WSCALE  AND  ZSCALE 
I F (ZSCALE. EQ.O. )  ZSCALE=1. 

I F (WSCALE . EQ. 0. )  WSCALE=1. 

C*»  READ  IN  HISTOGRAM  DISTANCE  AND  WEIGHT/LENGTH  POINTS  FROM  DATA2  FILE. 

C  DISTANCES  ARE  MEASURED  FROM  THE  FORWARD  PERPENDICULAR. 

120  READ (2,*)  IDIG,ZDIG,WXDIG 
I F ( EOF ( 2 ) . NE . 0 )  GO  TO  125 
LL=LL+1 

WA(1,LL)=ZDIG*ZSCALE 
WA ( 2, LL ) =WXD I G*WSCALE*COLW I D 
IWA(3,LL)=-1 
GO  TO  120 
125  NHI S=LL 
C 

C**  WSUM  AND  WMOM  ARE  CALCULATED  FROM  WEIGHT  COLUMNS  THAT  OVERLAP  SHIP 
C  SEGMENTS  OF  CROSS  SECTIONS  APPEARING  ON  DATA1 .  THIS  RESULTS  IN 

C  MISSING  SECTIONS  BEING  TREATED  AS  NEUTRAL.  HOWEVER,  EXTRA  WEIGHTS 

C  DUE  TO  FLOODING,  IF  PRESENT,  ARE  INCLUDED  BELOW  IN  WSUM  AND  WMOM. 

C  (  PROGRAM  ASSUMES  COLW ID  .LE.  SEGWID  ) 

D I SPL=WSUMX=WSUM=WMOM=0 . 

HLFCOL=COLW ID/2. 

HLFSEG=$EGW ID/2. 

DO  155  L=1 ,NH I S 
D I SPL=D I SPL+WA ( 2, L ) 

ZWM I N=WA ( 1 , L ) -HLFCOL 
ZWMAX=WA ( 1 , L ) +HLFCOL 
ICHK=0 

DO  150  1=1,111 
ZB=NS( I )*SEGWID 
ZBM I N=ZB-HLFSEG 
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ZBMAX=ZB+HLF  SEG 

I F ( ZWM I N . GE . ZBMAX . OR . ZWMAX . LE . ZBM I N )  GO  TO  150 
I F { ZWM I N . LT . ZBM I N )  GO  TO  130 
IF (ZWMAX. GT. ZBMAX)  GO  TO  140 
WSUM=WSUM+WA(2,L) 

WMOM=WMOM+WA ( 1 , L ) *WA ( 2 , L ) 

GO  TO  155 

1 30  WLUMP= ( ZWMAX-ZBM I N ) /COLW I D*WA ( 2, L ) 

ZLUMP= ( ZBM I N+ZWMAX )/2. 

GO  TO  145 

1 40  WLUMP= ( ZBMAX-ZWM I N ) /COLW I D*WA (2, L) 

ZLUMP=(ZWMIN+ZBMAX)/2. 

145  WSUM=WSUM+WLUMP 

WMOM=WMOM+WLUMP*ZLUMP 
IF( ICHK.EQ.1 )  GO  TO  155 
ICHK=1 

150  CONTINUE 
155  CONTINUE 

IF (.NOT. FLOODED)  GO  TO  175 
C 

C*  READ  IN  EXTRA  WEIGHT  AND  DISTANCE  POINTS 
READ  1 60,NXTRA 
160  FORMAT (110) 

DO  170  1=1, NXTRA 
LL=LL+1 

READ  1 65,ZXTRA, WXTRA 
165  FORMAT (2E1 0. 0) 

WA ( 1 , LL ) =ZXTRA 
WA(2,LL)=WXTRA 
WMOM=WMOM+WXTRA*ZXTRA 
WSUMX=WSUMX+WXTRA 
IWA(3,LL)=0 
170  CONTINUE 

WSUM=WSUM+WSUMX 
D I SPL=D I SPL+WSUMX 
C 

C*»*  BEGIN  ITERATION  TO  FIND  SLOPE  AND  YINT  FROM  WHICH  YCSO  ARRAY  IS  CALCULATED. 
C*  YCSO  GIVES  VERTICAL  DISTANCE  BETWEEN  CALCULATION  AND  SHIP  ORIGINS. 

C 

C*  ASSUME  AN  INITIAL  DRAFT  AND  SLOPE  FOR  LINE  CONNECTING  SECTION  WATERLINES. 

C  POWER  LAW  IS  FROM  F I T  TO  TYPICAL  SHIP  DRAFT-DISPLACEMENT  DATA. 

175  Y I  NT  = 1 .96*01 SPL**0. 26 
SLOPE=0. 

C*  MODIFY  YINT  VALUE  IF  UNITS  USED  APPEAR  TO  BE  METRIC. 

IF(DENSITY#35. .GT. 1 0. )  YINT=YINT/3. 

C 

180  BSUM=BMOM=BP=BPZ=BPZ2=0. 

DO  190  I S=1 , I  I  I 
I  I  =  I S 

ZB  I =FLOAT (NS( I  I ))*SEGWID 

YCSO( I  I )=Y I NT+SLOPE*ZB I -YDCO( I  I ) 

JJ  =  JJA( I  I ) 

KK=KKA( I  I ) 

CALL  AREACAL 
BI=AREA( I  I )*SEGDEN 
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BP  I =2.*(DCU( I  I )-XCS0( I  I ) )*SEGDEN 
BSUM=BSUM+B I 
BM0M=BM0M+B I *ZB I 
BP=BP+BPI 
BPZ=BPZ+BP I *ZB I 
BPZ2=BPZ2+BPI*ZBI**2 
WA( 1 ,  LL+I I )=ZB I 
WA(2,LL+I I )=BI 
IWA(3,LL+I I )  =  I  I 
190  CONTINUE 

D1=WSUM-BSUM 

D2=WM0M-BM0M 

IF (ABS(D1/WSUM)+ABS(D2/WM0M) .LT. 1 .E-9)  GO  TO  200 
DTERM= ( BP*BPZ2-BPZ**2 ) 

Y I NT=Y I NT+ ( D 1 *BPZ2-D2*BPZ ) /DTERM 
SLOPE=SLOPE+ ( D2*BP-D 1 *BPZ ) /DTERM 
GO  TO  180 
C 

C*»*  PRINT  EQUILIBRIUM  CALCULATION  RESULTS. 

200  AKEEL=ATAN(SL0PE)/PID1 80 
NWA=LL+ 1  I  I 

CALL  SORT (WA,NWA*3, 3) 

PRINT  210, (TITLE( J), J=1 ,8) 

210  FORMAT (*1*4X,8A10//5X*RESULTS  OF  SHIP  EQUILIBRIUM  CALCULATION 
1 (UNLISTED  STATIONS  CONSIDERED  NEUTRAL)*) 

PRINT  21 5,WSUMX 

215  FORMAT (*0*9X*EXTRA  WEIGHT  DUE  TO  FLOODING  =*F8.2) 

PRINT  220,01 SPL,AKEEL 

220  FORMAT (*0*9X*D I SPLACEMENT  =*F9.2//10X*ANGLE  OF  KEEL  =  *F5.2 
1*  DEGREES  FROM  HORIZONTAL  (POS  =  BOW  UP)*) 

PRINT  230 

230  FORMAT ( *0*// 1 4X*STAT I 0N*3X*FRAME*3X*FREEB0ARD*3X*DRAFT*3X 
1 *D I  ST  FROM  STA  0*3X*BU0YANCY*4X*W(DATA2)*) 

IF(NXTRA.GT.O)  PRINT  235 
235  FORMAT (*+*T95*W( EXTRA)*) 

DO  250  I =1 ,NWA 
K= I WA( 3, I ) 

IF (K.EQ.-1 )  PRINT  240, WA( 1,1) ,WA(2, I ) 

240  FORMAT (T57,F6.2,T82,F8.2) 

IF (K.EQ.O  )  PRINT  245, WA( 1 , ! ) ,WA(2, I ) 

245  FORMAT (T57.F6. 2,T95,F7. 2) 

IF(K.LE.O)  GO  TO  250 
JJA(K) 

FREEBD=RC IN(J,K)*SIN(PCIN(J,K)) -YCSO ( K ) 

DRAFT=YCSO ( K ) -YCL ( K ) 

PRINT  80, NS(K),FRAME(K),FREEBD, DRAFT, WA(1, I ),WA(2, I ) 

250  CONTINUE 
300  RETURN 
END 
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SUBROUTINE  SHI  PLOT (MODE) 

»*  SHI  PLOT  IS  A  GOULD  OR  CALCOMP  PLOTTING  SUBROUTINE.  MOST  OF  THE 
SUBROUTINE  CALLS  DEPEND  ON  A  PROPRIETARY  SOFTWARE  PACKAGE  FOR 
A  GOULD  ELECTROSTATIC  PLOTTER. 

SECTIONS  ARE  PLOTTED  ON  INDIVIDUAL  PLOTS  (MULTI =TRUE) 

OR  TOGETHER  ON  A  SINGLE  PLOT  (MULT I =FALSE ) . 

CONTENT  OF  EACH  PLOT  IS  DETERMINED  BY  FLAGS  SET  IN  THE  I  PLOT  ARRAY. 

I  PLOT ( I ) =0  ...SKIP  l-TH  PLOT  CATEGORY. 

I  PLOT ( 1 )=N  ...DRAW  FITTED  SPLINE  FUNCTION  CURVE.  (FOR  N  SEE  BELOW) 

I  PLOT (2)=N  ...PLOT  INPUT  POINTS  DEFINING  SECTION. 

I  PLOT ( 3 )=N  ...DRAW  FITTED  MAPPING  FUNCTION  CURVE. 

IPL0T(4)=N  ...PLOT  Z-PLANE  IMAGES  OF  POINTS  ON  UNIT  CIRCLE. 

N  IS  AN  INTEGER  FROM  1  TO  9  THAT  CONTROLS  THE  THICKNESS  OF  THE  PLOTTED 
POINT  OR  LINE.  A  VALUE  OF  4  PRODUCES  A  STANDARD  THICKNESS.  SMALLER 
VALUES  ARE  LIGHTER  AND  LARGER  VALUES  ARE  DARKER. 

FMAG  IS  A  MAGNIFICATION  FACTOR  THAT  AFFECTS  ALL  PLOT  DIMENSIONS  (1.-1  TO  1) 
XEND  AND  YEND  ARE  LENGTHS  OF  THE  X  AND  Y  AXES  RELATIVE  TO  THE  SHIP  ORIGIN. 
THESE  MUST  BE  GIVEN  IN  TENS  OF  FEET  (E.G. , 10,20,30,40) . 

**  ALL  SHIPFIT  CALLS  TO  PLOTTING  SUBROUTINES  ARE  CONTAINED  IN  SHI  PLOT. 

**  IT  IS  CALLED  BY  SHIPFIT  WITH  THE  FOLLOWING  ARGUMENT  VALUES, 

M0DE=1  ...  INITIALIZE  PLOTS. 

MODE =2  ...PLOT  I l-TH  SECTION  (ACCORDING  TO  IPLOT  VARIABLE). 

M0DE=3  ...END  ALL  PLOTTING. 

COMMON  /BL0CK1 /  SUBMERG,NS(24 ) ,FRAME(24 ) ,KKA(24 ) ,KK, J JA(24 ) , J J, 

1  XCSO ( 24 ) , YCSO ( 24 ) , RC I N ( 64 , 24 ) , PC  I N ( 64 , 24 ) , CF ( 28 , 24 ) , E ( 24 ) , DCU ( 24 ) 
Z,YCL(24) , AREA (24) .DENSITY, P,R,RP,RPP, 11,111 
COMMON  /BL0CK2/  AA( 26, 3) ,TI TLE(8) ,XSHI FT, YSHI FT, ALPHA, BETA, 

1  CX( 24,64 ) ,CY( 24,64 ), MM, NN.NSEGS 

COMMON  /BL0CK3/  P I D1 80, YDC0(24 ) ,CA( 28, 24 ) ,RC2( 64, 24 ) ,X(9, 64 ) 

DIMENSION  XPLT ( 200 ) , YPLT ( 200 ) , I  PLOT ( 1 0 ) , P I N ( 9 ) 

LOGICAL  SUBMERG, MULTI 

DATA  PIN/. 125,. 25,. 5,1., 2., 3., 4., 6., 9./ 

DATA  BCDX,BCDY/1 HX, 1 HY/ 

DATA  P 1 02/ 1 . 5707963267949/ 

I F (MODE. EQ. 1 )  GO  TO  1 
IF (M0DE.EQ.2)  GO  TO  35 

**  END  ALL  PLOTTING. 

CALL  PEND 
GO  TO  100 

*»  INITIALIZE  PLOT  PARAMETERS  AND  DRAW  AXES  IF  SINGLE  PLOT. 

READ  (AND  SET  DEFAULT  VALUES  FOR)  FMAG,  XEND,  AND  YEND. 

1  READ  5, (IPLOT(I), 1=1, 10), MULTI, FMAG, XEND, YEND 
5  FORMAT (1011,110, 3E 10.0) 

IF(FMAG.EQ.O. )  FMAG= 1 . 

IF(XEND.EQ.O. )  XEND=30. 
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IF (YEND.EQ.O. )  YEND=30. 

C 

NXD I V=XEND/ 1 0. +. 5 
NYDI V=YEN0/10.+. 5 
NPLT=1 98 

XPLT ( 1 99 )=YPLT ( 1 99 )=0. 

THSTP=PID2/FL0AT(NPLT-1 ) 

IF (MULT I )  GO  TO  20 
C 

C**  SET  PARAMETERS  FOR  SINGLE  PLOT  AND  DRAW  AXES. 

CALL  PS  I ZE ( -1 .,-0.5,16., 10., -1 ) 

CALL  FACTOR ( FMAG*2 . , FMAG*2. ) 

I FL I P= I  I  1/2 
IF (SUBMERG)  GO  TO  10 

C*  SET  LENGTH  OF  NEGATIVE  Y  AXIS  (REL.TO  SO)  FOR  SURFACE  SHIP. 
YL=6. 

IS0=1 

Y0Y=2. 

GO  TO  15 

C*  SET  LENGTH  OF  NEGATIVE  Y  AXIS  (REL.TO  SO)  FOR  SUBMERSIBLE. 

10  YL=3. 

I  S0=2 
Y0Y=2. 

THSTP=2. *THSTP 
15  YOX=YOY+YL* I  SO 
TENFT=YL/NYD I  V 
XL=NXD I V*TENFT 
XPLT ( 200 )=XEND/XL 
YPLT ( 200 ) =YEND/ YL 
C 

C*  DRAW  X  AND  Y  AXES  FOR  SINGLE  PLOT. 

CALL  I NTENS ( 2. ) 

CALL  MOVEA(0. , YOX) 

CALL  ANGLE (0.) 

CALL  AXCCW 

CALL  AXDEF ( 2*NXD I V*TENFT, 2*NXD I  V ) 

CALL  TXTS I Z( . 1 , . 05, . 1 ) 

CALL  TXTFMT ( °F°, 3,0) 

CALL  AXVAL ( -XEND, 1 0. ) 

CALL  MOVEA (XL, YOY ) 

CALL  YAX I S ( YL* I  SO, I SO*NYD I  V ) 

CALL  SYMBOL ( XL-3. , YOY-. 5,. 14, TITLE, 0., 80) 

CALL  REORGA(XL,YOY+YL) 

GO  TO  100 
C 

C*»  SET  PARAMETERS  FOR  MULTIPLE  PLOTS. 

20  IF (SUBMERG)  GO  TO  25 
C 

C*  SET  LENGTH  OF  NEGATIVE  Y  AXIS  (REL.TO  SO)  FOR  SURFACE  SHIP. 
YL=6. 

I  S0=  1 
Y0Y=2. 

GO  TO  30 
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C 

C*  SET  LENGTH  OF  NEGATIVE  Y  AXIS  (REL.TO  SO)  FOR  SUBMERSIBLE 
25  YL=3. 

I  S0=2 
Y0Y=2. 

THSTP=2. *THSTP 
30  YOX=YOY+YL* I  SO 
TENFT=YL/NYDI V 
XL=TEKFT*NXDIV 
XPLT ( 200 ) =XEND/XL 
YPLT ( 200 ) =YEND/ YL 
GO  TO  100 


C*»*  CREATE  AND  PLOT  DATA  ARRAYS  (DRAW  AXES  IF  MULTIPLE  PLOTS) 
35  IF(. NOT. MULTI )  GO  TO  40 
C 

C*»  DRAW  AND  LABEL  AXES  FOR  ll-TH  SECTION. 

CALL  PS  I ZEC-1 . ,-0. 5, 1 6. , 1 0. ,-1 ) 

CALL  F ACTOR (FMAG*2.,FMAG*2. ) 

CALL  INTENS(2. ) 

CALL  MOVEAtO. ,  YOX) 

CALL  ANGLE (0.) 

CALL  AXCCW 

CALL  AXDEF (XL,NXD I  V ) 

CALL  TXTSIZ( . 1 , .05, .  1  ) 

CALL  TXTFMT ( °F°, 3, 0) 

CALL  AXVAL (0. ,10.) 

CALL  TXTSI Z( .  1 5, . 1 , . 1 5) 

CALL  AXLAB( 1 ,BCDX) 

CALL  MOVEACO. , YOY) 

CALL  YAX I S ( YL* I  SO, I SO*NYD I  V ) 

CALL  TXTS I Z( . 1 , . 05, . 1 ) 

CALL  TXTFMT (°F°, 3,0) 

CALL  AXVAL ( -YEND, 1 0. ) 

CALL  TXTS IZ( . 1 5, . 1 , . 1 5) 

CALL  AXLAB(1,BCDY) 

CALL  SYMBOL(0. , YOY-. 5,  . 14, TITLE, 0. ,70) 

CALL  SYMBOL (XL-2. 5, YOY+. 5, . 1 4, 1 4HSTAT ION  NUMBER, 0. , 1 4 ) 
CALL  NUMBER(XL-. 4, YOY+. 5, . 1 4, FLOAT (NS ( I  I )),0.,0) 

CALL  REORGA(0. ,YOY+YL) 

C 

C»*  PLOT  DATA  FOR  ll-TH  SECTION. 

40  XS I GN= 1 . 

IF(. NOT. MULTI .AND. I  I .GT. IFLIP)  XSIGN=-1 . 

I F ( (PLOT ( 1 ) . EQ.O)  GO  TO  60 
C 

C*  DRAW  SPLINE  FUNCTION  CURVE  (FITTED  TO  INPUT  POINTS). 

I F ( SUBMERG )  GO  TO  45 
PL  =ATAN2 ( YCL (II), XCSO (II)) 

PU=ATAN2(YCS0(  I  I  ) ,  DCU  ( ID) 

GO  TO  50 

45  PL=ATAN2(YCL(  i  I  ),XCSO(  ID) 

PU=ATAN2(DCU(  I  D,XCSO(  ID) 

50  PRANGE=PU-PL 
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PSTP=PRANGE/FL0AT(NPLT-1 ) 

P=PL-PSTP 
DO  55  I =1 ,  NPLT 
P=P+PS TP 
CALL  SPLCAL 
C=COS(P) 

S=S I GN ( SQRT ( 1 . -C*C ) ,  P ) 

XPLT ( I )=(R*C-XCSO( I  I ) )*XS IGN 
YPLT ( I )=(R*S-YCSO( II)) 

55  CONTINUE 

CALL  I NTENS ( P I N ( I  PLOT ( 1 ) ) ) 

CALL  L I NE (XPLT, YPLT, NPLT, 1,0,0) 

60  IF ( I  PLOT (2) .EQ.O)  GO  TO  75 

PLOT  INPUT  POINTS  DEFINING  ll-TH  SECTION. 

DO  65  J=1,JJ 
NJS=J 

C=COS(PCIN(  J,  ID) 

S=S  I  GN  ( SQRT  ( 1 .  -C*C ) ,  PC  I N  ( J ,  ID) 

XPLT ( J )  =  (RC I N( J, I  I )*C-XCSO( I l))*XSIGN 
YPLT ( J)  =  (RCIN( J, I  I )*S-YCSO< I D) 

IF (SUBMERG)  GO  TO  65 
IF (YPLT ( J ) . LE. 0. )  GO  TO  65 
NJS= J-1 
GO  TO  70 
65  CONTINUE 

70  XPLT(NJS+1 )=YPLT (NJS+1 )=0. 

XPLT (NJS+2)=XEND/XL 
YPLT (NJS+2)=YEND/YL 
CALL  I NTENS ( P I N ( I  PLOT ( 2 ) ) ) 

CALL  L I NE (XPLT, YPLT, NJS, 1 , ~1 , 1 ) 

75  IF ( I  PLOT ( 3 ) .EQ.O)  GO  TO  92 

DRAW  MAPPING  FUNCTION  CURVE  FOR  ll-TH  SECTION. 

DO  80  M=1 ,NPLT 
XPLT (M)=0. 

YPLT(M)=0. 

80  CONTINUE 

DO  90  M=1 ,NPLT 

THETAM=THSTP»FL0AT(M-1 )-PID2 
DO  90  N=1,NN 

ARG= ( ALPHA -BE  TA  *N ) *THE TAM 

IF ( SUBMERG. AND. MOD (N, 2) .EQ.O)  GO  TO  85 

XPLT (M)=XPLT (M)+AA(N, 1 )*COS(ARG) 

YPLT(M)=YPLT(M)+AA(N,  1 )*S IN( ARG) 

GO  TO  90 

85  XPLT (M)=XPLT (M)-AA(N, 1 )*SIN(ARG> 

YPLT (M)=YPLT(M)+AA(N, 1 )*COS ( ARG) 

90  CONTINUE 

CALL  I NTENS ( P I N ( I  PLOT ( 3 ) ) ) 

CALL  LINE (XPLT, YPLT, NPLT, 1,0,0) 

92  IF ( I  PLOT (4) .EQ.O)  GO  TO  100 

PLOT  MAPPING  FUNCTION  IMAGE  POINTS  FOR  THE  ll-TH  SECTION 
DO  95  M=1,MM 
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XPLT (M)=X(6,M) 

YPLT (M)=X(7,M) 

95  CONTINUE 

XPLT (MM+1 )=YPLT (MM+1 )=0. 

XPLT ( MM+2 ) =XEND/XL 
YPLT ( MM+2 ) =YEND/ YL 
CALL  I NTENS ( P I N ( I  PLOT ( 4 ) ) ) 

CALL  LINE (XPLT, YPLT, MM, 1,-1, 2) 
100  RETURN 
END 
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SUBROUTINE  SORT (A,NA,LG) 

THIS  IS  A  SIMPLE  INSERTION  SORT 

DIMENSION  A(LG,1 ),T(25) 

C 

NG=NA/LG 

I F{NG  .LT.  2)  RETURN 
C 

DO  150  I  =  2,NG 

l F ( A ( 1 , I ) ,GE.A( 1 , 1-1 ) )  GO  TO  150 
J=  I -1 

50  DO  60  K=1 , LG 
T(K)=A(K, I ) 

A(K, J+1 )=A(K, J) 

60  CONTINUE 
J=J-1 

C  THE  TEST  FOR  (J  .EQ.  0)  IS  FASTER  THAN  A  TEST  FOR  (J  .EQ.  1) 
IF ( J  .EQ.  0)  GO  TO  100 
IF(T(1 )  .LT.  A( 1 , J ) )  GO  TO  50 
C  A  SUBSCRIPT  OF  THE  FORM  (0+1)  IS  ACCEPTABLE. 

100  DO  110  K=1,LG 
A(K, J+1 )=T(K) 

110  CONTINUE 
150  CONTINUE 
RETURN 
END 
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SUBROUTINE  SPLCAL 

COMMON  /BLOCK 1 /  SUBMERG,NS(24),FRAME(24),KKA(24),KK, JJA(24), JJ, 

1  XCSO ( 24 ) ,  YCSO  ( 24 ) ,  RC I N  ( 64 , 24 ) ,  PC  I N  ( 64 , 24 ) ,  CF  ( 28 , 24 ) ,  E  ( 24 ) , DCU ( 24 ) 
2, YCL ( 24 ) , AREA ( 24 ) , DENS  I TY, P, R, RP, RPP, 11,111 
I J=KK-1 

SUM1 =SUM2=SUM3=0. 

DO  2  K=1 , J J,  I  J 

I F (P.LE.PC I N(K, I  I ))  GO  TO  3 

CONTINUE 

K=JJ 

IF(K.LE.KK)  GO  TO  5 

K=(K-KK)/I J 

DO  4  J=1 ,K 

DP=PCIN<J*I J+1,  I  I )-P 

C341 =(CF( J+3, I  I )-CF ( J+4, I  I ))*DP 

C342=C341»DP 

C343=C342*DP 

SUM1=SUM1+C341 

SUM2=SUM2+C342 

SUM3=SUM3+C343 

RPP=(2.*CF(3, 1 1 ) )+6.*( (P*CF(4, 1 1 ))+SUM1 ) 

RP=CF (2, I  I )+P*( (2.*CF (3, 1  I ) )+3.*(P*CF(4, I  I ) ) )-3.*SUM2 
R=CF (1,1  I )+P»(CF(2, 1 1 )+P*(CF (3, 1  I )  +  (P*CF(4, 1  I ))))+SUM3 
RETURN 
END 
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SUBROUTINE  SPLINT  <N,X,Y, JJ,A,R,S,C) 

DIMENSION  X(1),Y(1),C(25) 

I J=JJ-1 
13  SUM=0. 

DO  500  K=1 ,N, I J 

IF  ( S-X(K) )  400,400,500 

500  CONTINUE 
K=N 

400  IF(K-JJ)  200,200,600 

600  K=(K-J J )/l J 

300  DO  800  J = 1 ,  K 
!K=J* I J+1 

800  SUM=SUM+(C(J+3)-C(J+4))*(X( IK)-S)**4 

200  SUM1 =0. 

DO  501  K=1 , N,  I  J 
IF  (R-X(K) )  401,401,501 

501  CONTINUE 
K=N 

401  IF(K-JJ)  201,201,601 

601  K=(K-JJ)/I J 

301  DO  801  J=1,K 
IK=J* I J+1 

801  SUM1=SUM1+(C(J+3)-C(J+4))*(X(IK)-R)**4 

201  A=(C( 1 )+C(2)*S/2. +C(3 ) *S**2/3.+C(4 )*S**3/4. )#S-(C( 1 )+C(2 )*R/2. 
1+C(3)*R**2/3.+C(4)*R^3/4. ) *R- . 25* ( SUM-SUM  1  ) 

11  RETURN 
END 
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SUBROUTINE  SPLSQ1  (NPT,X,Y,C,E,  JJ) 

C 

C#*»  SPLSQ1  FITS  A  CUBIC  SPLINE  FUNCTION  TO  INPUT  DATA  BY  A  LEAST  SQUARES 
C  TECHNIQUE.  IT  WAS  WRITTEN  BY  CHARLES  WEBER  AND  IS  AVAILABLE  THROUGH 
C  THE  NSWC/WHITE  OAK  PROGRAM  LIBRARY. 

C 

DIMENSION  DF(25),X(1 ),Y(1 ),C(1),V(25,25) 

I J= J J-1 

NR=(NPT-1 )/l J+3 
NUFSED=0 
DO  1  1=1, NR 
C(l)=0. 

DO  1  J=1 ,NR 

1  V(l,J)=0. 

8  DO  10  1=1, NR 
10  DF ( I )=0. 

DO  20  1=1, NPT 
DF ( 1  )  =  1 . 

DF (2)=X(  I ) 

DF ( 3 ) =X ( I )**2 
DF (4)=DF(3)*X( I ) 

IF(I.LE.JJ)  GO  TO  6 

4  DF(4)=DF (4)+(X(JJ)-X(l ) )**3 
N= I  —  ( I *( JJ-2J+1 )/ 1 J 

IF (N.GT.NR-3)  N=NR-3 

IK= I J*(N-1 )+1 

DF ( N+3 ) =- ( X ( I K ) -X ( I ) ) #*3 

IF (N.EQ.2)  GO  TO  6 

N2=N+2 

DO  5  K=5,N2 

IK=(K-4)*I J+1 

IL=IK+I J 

5  DF(K)=-(X( IK)-X( I ) )**3+<X( IL)-X< l))»*3 

6  IF (NUFSED.EQ.O)  GO  TO  7 

Y1  =0. 

DO  3  J=1,NR 
3  Y1 =Y1 +C( J )#DF ( J ) 

SUM=SUM+(Y1-Y( I ))**2 
GO  TO  20 

7  DO  2  J=1 ,NR 

C ( J ) =C ( J ) +Y ( I ) »DF ( J ) 

DO  2  K=1 ,NR 

2  V(K,J)=V(K,J)+DF(K)*DF(J) 

20  CONTINUE 

IF(NUFSED.EQ. 1 )  GO  TO  9 
CALL  MATRIX  ( V,NR, 25, C, 1 ,DE, ID) 

SUM=0. 

NUFSED=1 

IF(E.LT.O. )  RETURN 
GO  TO  8 

9  E=SQRT( SUM/FLOAT (NPT-NR ) ) 

RETURN 

END 
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SUBROUTINE  SUBBEAM (B02) 

C**»  SU8BEAM  CALCULATES  THE  HALF  BEAM  OF  A  SUBMERSIBLE  FORM 

COMMON  /BL0CK1 /  SUBMERG, NS ( 24 ) , FRAME ( 24 ) , KKA ( 24 ) ,  KK,  J  JA ( 24 ) ,  J  J , 

1  XCSO ( 24 ) , YCSO ( 24 ) , RC I N ( 64 , 24 ) , PC  I N ( 64 , 24 ) , CF ( 28, 24 ) , E ( 24 ) ,  DCU ( 24 ) 
2, YCL(24),AREA(24), DENSITY, P,R,RP,RPP, 11,111 
XCMAX=0. 

DO  10  J=1,JJ 

XC=RC  I N  ( J ,  I  I  )*COS(PCIN(J,  ID) 

XCMAX=AMAX1 (XCMAX.XC) 

IF (XCMAX.EQ.XC)  JMAX=J 
10  CONTINUE 

P=PCIN( JMAX, I  I ) 

20  CALL  SPLCAL 
C=COS(P) 

T=TAN(P) 

F=RP-R*T 

IF (F. EQ.O. )  GO  TO  30 
DELP=F/(RPP-RP*T-R/C**2) 

P=P-DELP 

IF (DELP/P.LT. 1 .E-6)  30,20 
30  BD2=R*C-XCS0( I  I ) 

RETURN 

END 
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SUBROUTINE  XLIMIT(PIN.XC) 

C**»  XL  I  MIT  FINDS  BY  ITERATION  THE  INTERSECTION  OF  THE  SHIP  FORM  CURVE  AND  THE 
C  ABSCISSA  OF  THE  SHIP  COORDINATE  SYSTEM.  OUTPUT  IS  RELATIVE  TO  CO. 

COMMON  /BLOCK!/  SUBMERG,NS(24),FRAME(24),KKA(24),KK, JJA(24), JJ, 

1  XCSO ( 24 ) , YCSO { 24 ) , RC I N ( 64 , 24 ) , PC  I N ( 64 , 24 ) , CF ( 28 , 24 ) , E ( 24 ) , DCU ( 24 ) 

2, YCL(24),AREA(24), DENSITY, P,R,RP,RPP, 11,111 
P=PIN 

10  CALL  SPLCAL 
S=SIN(P) 

C=SQRT(1.-S*S) 

F=R*S-YCSO( I  I ) 

I F  ( F .  EQ .  0 . )  GO  TO  20 
DELP=F/(RP*S+R*C) 

P=P-DELP 

I F ( ABS ( DELP/P ) . LT . 1 . E-6 )  20,10 
20  XC=R*C 
RETURN 
END 


SUBROUTINE  YLIMIT(PIN.YC) 

C***  YLIMIT  FINDS  BY  ITERATION  THE  INTERSECTION  OF  THE  SHIP  FORM  CURVE  AND  THE 
C  ORDINATE  OF  THE  SHIP  COORDINATE  SYSTEM  IN  THE  VICINITY  OF  THE  CO  POINT 
C  (PI N,R(PI N ) ) .  OUTPUT  IS  RELATIVE  TO  CO. 

COMMON  /BL0CK1/  SU8MERG, NS(24) ,FRAME(24) ,KKA(24) ,KK, JJAC24) , J J, 

1  XCSO ( 24 ) , YCSO ( 24 ) , RC I N ( 64 , 24 ) , PC  I N ( 64 , 24 ) , CF ( 28 , 24 ) , E ( 24 ) , DCU ( 24 ) 

2, YCL ( 24 ) , AREA ( 24 ) , DENS  I TY, P, R,RP,RPP, 11,111 
P=PI  N 

10  CALL  SPLCAL 
C=COS(P) 

$=S I GN ( SORT ( 1 . -C*C ) , P ) 

F=R*C-XCSO( 1 1 ) 

IF (F.EQ.O. )  GO  TO  20 
DELP=F/ (RP*C-R*S ) 

P=P-DELP 

I F ( ABS ( DELP/P ) . LT . 1 . E-6 )  20,10 
20  YC=R*S 
RETURN 
END 
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